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Finite  Element  Modeling  of  Tire-Terrain  Interaction 


SALLY  A.  SHOOP 


INTRODUCTION 

Vehicle  mobility  on  unpaved  surfaces  is  important 
to  the  military  as  well  as  to  agriculture,  forestry,  min¬ 
ing,  and  construction  industries.  The  problems  can  be 
grouped  into  two  major  categories:  predicting  vehicle 
performance  on  various  terrains  (will  it  get  stuck; 
how  much  traction  or  pull  is  available  to  climb  or 
pull)  and  estimating  the  consequences  of  the  vehicle 
passage  (rut  formation,  shearing/tearing  of  roots,  soil 
compaction,  and  the  effects  of  these  on  vegetation 
and  erosion). 

The  objective  of  this  project  is  to  produce  a  three- 
dimensional  finite  element  model  of  tire-terrain  inter¬ 
action  that  can  be  used  to  explore  the  effects  of  tire 
and  terrain  variables  on  vehicle  mobility  (traction  and 
motion  resistance)  and  terrain  deformation  (rut  for¬ 
mation  and  shearing).  Such  a  model  would  enable 
detailed  analysis  of  the  complex  interactions  resulting 
from  contact  friction  and  would  further  the  under¬ 
standing  of  off-road  vehicle  mobility  by  defining 
critical  mechanisms  involved  in  vehicle  traction  and 
motion  resistance.  Ultimately  the  model  generated 
would  be  used  for  tire  design  and  specification  for 
off-road  vehicles  (for  construction,  mining,  and  rec¬ 
reation),  for  vehicle  performance  prediction,  and  for 
terrain  damage  prediction  and  reduction  of  the 
environmental  impact  of  off-road  travel.  Previously  a 
three-dimensional  simulation  of  contact  between  a 
deformable  tire  and  deformable  terrain  had  been  too 
difficult  and  computationally  time  consuming.  Re¬ 
cent  advancements  in  the  contact  formulations  of 
general-purpose  finite  element  codes  (e.g.  ABAQUS, 
HKS  1998)  and  increases  in  computer  processing 
speeds  have  brought  such  a  model  into  the  realm  of 
possibility. 

The  originality  of  this  research  lies  in  accounting 
for  the  deformable  nature  of  both  the  tire  and  the  ter¬ 
rain  in  a  fully  three-dimensional  model.  Both  tire  and 


terrain  are  highly  deformable.  A  rigid  approximation 
of  either  is  an  oversimplification  of  the  problem  and 
would  yield  erroneous  contact  conditions,  which  are 
critical  to  obtaining  correct  solutions.  The  incorpora¬ 
tion  of  the  third  (lateral)  dimension  is  similarly  nec¬ 
essary,  as  the  lateral  movement  of  the  terrain  (in  ad¬ 
dition  to  deformation  in  front  and  beneath  the  tire)  is 
readily  apparent  for  most  soft  terrain  materials  (soils 
and  snow).  The  effect  of  the  lateral  dimension  is  also 
evident  in  equations  for  rolling  resistance,  which  is  a 
function  of  the  width  of  the  rolling  object  (Bowden 
and  Tabor  1964,  Ludema  1996a).  Similarly,  empiri¬ 
cal  sinkage  and  motion  resistance  equations  by  Bek- 
ker  (1969)  for  soil  and  Richmond  (1995)  for  snow 
include  the  tire  or  track  width,  indicating  the  impact 
of  the  third  dimension  on  the  performance  of  the  tire 
or  track. 

The  details  of  the  tire-soil  modeling  problem  are 
divided  into  three  topic  areas:  1)  material  models  for 
the  terrain  material,  2)  tire  models  for  use  on  a  de¬ 
formable  substrate,  and  3)  the  combined  tire-terrain 
model  and  the  treatment  of  the  interface.  These  are 
presented  after  a  brief  historical  perspective  on  tire- 
terrain  interaction. 


BACKGROUND 

Off-road  vehicle  performance  prediction 

Off-road  vehicle  performance  analysis  is  a  three- 
dimensional,  nonlinear,  dynamic  problem.  Major 
early  efforts  were  led  by  Bekker,  at  the  University  of 
Michigan  and  the  U.S.  Army  Land  Locomotion 
Laboratory  (Bekker  1956,  1960,  1969),  and  his  con¬ 
temporaries  at  the  National  Tillage  Laboratory  at 
Auburn,  Alabama  (Gill  and  Vanden  Berg  1967). 
Karafiath  and  Nowatzki  (1978)  extended  the  state  of 
the  art  to  include  modeling  of  soil  deformation  under 
the  wheel  or  track  based  on  theoretical  soil  mechan- 
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ics  using  plasticity  theory.  They  developed  a  two- 
dimensional  finite  difference  model  of  tire-terrain 
interaction  that  simulated  the  plastic  soil  deformation 
under  a  tire,  accurately  reproducing  the  experimental 
measurements  of  sand  deformation  under  a  rigid 
wheel  by  Wong  and  Reece  (1967),  as  shown  in  Fig¬ 
ure  1.  More-recent  seminal  books  on  the  subject  are 
those  by  Yong  et  al.  (1984),  which  includes  some  of 
the  early  finite  element  modeling  of  tire-terrain  inter¬ 
action,  and  Wong  (1989),  who  approaches  numerical 
modeling  of  tracked  vehicles  based  on  Bekker’s 
semi-empirical  representation  of  the  terrain. 


Figure  1.  Sand  deformation  under  a  towed  wheei 
moving  to  the  right.  Dashed  iines  show  the  soii 
faiiure  zones.  (From  Wong  and  Reece  1967.) 


Because  vehicle-terrain  interaction  is  compli¬ 
cated,  comprehensive  modeling  of  off-road  mobility 
is  usually  done  using  empirical  algorithms  based  on 
large  data  sets,  such  as  the  NATO  Reference  Mobil¬ 
ity  Model  (Ahlvin  and  Haley  1992).  These  are  sup¬ 
plemented  with  high-fidelity  simulations  using  so¬ 
phisticated  vehicle  dynamics  models,  but  in  the 
future  these  simulations  will  use  vehicle-terrain  nu¬ 
merical  models. 

Vehicle  movement  on  snow 

Vehicle  mobility  on  snow  has  some  unique  char¬ 
acteristics  and  terminology.  The  performance  of  a 
vehicle  in  snow  is  usually  limited  by  the  bearing 
strength  of  the  snow,  which  controls  the  amount  of 
sinkage.  Sinkage  is  a  function  of  the  compressibility 
and  shear  strength  of  the  snow  in  response  to  the  ve¬ 
hicle  load.  Since  snow  strength  is  difficult  to  measure 
in  the  field  on  a  routine  basis,  vehicle  performance  is 
often  related  to  index  parameters  that  are  easier  to 
measure,  such  as  the  snow  depth  and  density. 


In  this  study,  only  fresh  snow  within  a  limited 
range  of  initial  density  (less  than  250  kg/m^)  was 
modeled.  The  snow  depth,  however,  ranged  from 
“shallow”  to  “deep.”  Shallow  snow  is  defined  as 
snow  where  the  “pressure  bulb”  (the  snow  compacted 
due  to  the  weight  of  the  vehicle)  intersects  a  rigid 
interface.  In  deep  snow  the  pressure  bulb  does  not 
intersect  a  rigid  interface  below  the  snow  (Fig.  2). 


Figure  2.  Shallow  and  deep  snow  under  a  wheeled 
vehicle. 


Since  the  pressure  bulb  forms  in  response  to  a  load 
on  the  snow,  the  classification  of  deep  or  shallow 
snow  is  vehicle  dependent  and  therefore  is  a  function 
of  the  vehicle’s  weight,  load  distribution,  tire  or  track 
geometry,  tire  pressure,  etc.  Nonetheless,  a  rule  of 
thumb  is  that  snow  deeper  than  two  thirds  of  the 
wheel  radius  is  usually  “deep”  snow.  Except  for  very 
specialized  vehicles,  greater  snow  depths  will  impede 
wheeled  vehicle  travel  because  of  excessive  sinkage. 
Generally,  tracked  vehicles  with  low  ground  contact 
pressure  are  used  to  traverse  deep  snow.  Figure  3 
shows  cross  sections  of  the  pressure  bulb  under  a 
vehicle  for  shallow  and  deep  snow. 

TIre-terraIn  modeling 

Because  of  the  complexities  involved  in  vehicle- 
terrain  interaction  and  the  desire  to  incorporate  more 
of  the  physics  into  performance  models,  there  is  great 
interest  in  applying  numerical  modeling  techniques. 
Several  researchers  have  produced  numerical  models 
simulating  tire-terrain  interaction  over  the  last  20 
years.  Many  of  the  models  are  two-dimensional 
simulations  of  a  rigid  wheel  on  a  deformable  surface. 
Yong  and  Fattah  (1976)  and  Yong  et  al.  (1978)  intro¬ 
duced  the  first  finite  element  models  of  tire-terrain 
interaction  using  a  two-dimensional  rigid  tire  on  de¬ 
formable  terrain  and  then  progressing  to  an  elastic 
tire  and  calculating  the  contact  using  Hertz  contact 
theory.  This  model  was  difficult  to  use,  as  it  needed  the 
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a.  Snow  deformation  under  a  Jeep  Cherokee  (2.7 
ton,  24  kN)  in  19-cm  “shaiiow”  snow. 


b.  Snow  deformation  under  a  BV202  (a  3.5-ton,  31- 
kN,  articuiated  tracked  vehicie)  in  a  122-cm 
“deep”  snow  cover.  (From  Harrison  1975.) 


Figure  3.  Snow  cross  sections  showing  pressure  buib  formation. 


contact  displacements  or  stresses  as  input  values. 
Nevertheless,  it  was  a  major  advance  in  applying 
numerical  methods  to  tire-terrain  interaction.  Model¬ 
ing  the  interface  between  a  deformable  soil  and  a 
deformable  wheel  has  been  both  problematic  and 
computationally  time  consuming.  Pi  (1988)  modeled 
a  two-dimensional  elastic  wheel  (connected  springs) 
on  a  viscoelastic  soil  for  high-speed  landing  of  air¬ 
craft  on  soil.  Others  have  used  estimates  of  the  con¬ 
tact  stress  distribution  predefined  on  the  soil  surface 
(Saliba  1990,  Chi  and  Tessier  1995).  This  approach 
simplifies  the  problem  and  is  suitable  for  analyzing 
different  terrain  parameters  but  has  the  major  disad¬ 
vantage  of  requiring  an  estimate  of  the  contact  dis¬ 
placements  or  stresses  a  priori  when  this  is  nearly 
impossible  to  predict  and  very  difficult  to  measure. 

Continued  advances  in  computing  capabilities  and 
in  general-purpose  finite  element  codes  have  enabled 
researchers  to  concentrate  on  the  physics  of  the 
model  rather  than  on  code  development.  Sophisti¬ 
cated  two-dimensional  models  of  a  rigid  wheel  on 
deformable  soil  have  been  described  by  Liu  and 
Wong  (1996)  and  Foster  et  al.  (1995).  Models  have 
recently  been  extended  to  three  dimensions  by  Chi,* 
Chiroux  et  al.  (1997),  and  Shoop  et  al.  (1999). 

To  date,  the  most  realistic  representation  of  a 
pneumatic  tire  on  deformable  soil  was  done  in  two 
dimensions  at  IKK  (Automotive  Research  Institute) 
in  Hamburg,  Germany  (Aubel  1993,  1994)  (Fig.  4). 
Fervers  (1994)  extended  AubeTs  model  to  study  the 
effect  of  lug  design  on  tire  performance  (Fig.  5).  A 
good  review  of  the  tire-terrain  research  program  at 


IKK,  including  the  developments  by  Aubel  and 
Fervers,  is  given  in  Schmid  (1995).  Fervers  further 
developed  his  model  to  a  pseudo-three-dimensional 
representation  and  used  it  in  several  applications, 
including  travel  over  rough  surfaces,  soil  compaction 
beneath  a  rolling  tire,  and  the  influence  of  tread  de¬ 
sign  and  slip  on  tire  performance  (Fervers  1997, 
1999a,  1999b). 

Current  efforts  toward  the  development  of  a  three- 
dimensional  model  of  a  deformable  tire  on  deform¬ 
able  terrain  are  being  undertaken  by  1)  the  U.S.  Army 
Cold  Regions  Research  and  Engineering  Laboratory 
(CRREL),  Goodyear,  and  Caterpillar  through  a  Co¬ 
operative  Research  and  Development  Agreement;  2) 
Tordesillas  (1996),  using  a  mathematical  contact  me¬ 
chanics  approach  at  the  University  of  Melbourne;  and 
3)  Shoop,  using  the  Darnell  tire  model  at  the  Univer¬ 
sity  of  Michigan,  Automotive  Research  Center,  and 
CRREL  (Mousseau  and  Hulbert  1996,  Darnell  et  al. 
1997,  Alverez  Sanz  1999,  Darnell,  in  progress).  The 
Automotive  Research  Center  (IKK)  in  Hamburg, 
Germany,  and  the  Transport  Technology  Research 
Laboratory,  Carleton  University,  Ottawa,  Canada, 
may  also  be  continuing  their  research.  Researchers  at 
the  Virtual  Proving  Ground  at  the  National  Automo¬ 
tive  Dynamic  Simulator  in  Iowa  have  recently  begun 
to  pursue  this  area  as  well.  The  research  described  in 
this  report  uses  the  Darnell  tire  model,  which  is  effi¬ 
cient  enough  to  make  simulations  of  a  rolling  tire 
achievable  in  near-real  time,  so  the  addition  of  a  de¬ 
formable  terrain  is  computationally  feasible.  This 
study  also  includes  the  development  of  a  terrain  ma¬ 
terial  model  for  snow. 
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*  Personal  communication,  L.  Chi,  Caterpillar,  Inc., 
Peoria,  IL,  1996. 


Figure  5.  Modeling  of  the  effect  of  lug  design  on  tread-terrain 
interaction  in  two  dimensions  by  Fervers  (1994). 


Tire  terminology 

Because  of  the  interdisciplinary  nature  of  this 
work,  some  of  the  tire  terminology  and  conventions 
are  presented  here,  as  they  are  used  throughout  the 
text.  The  tire  direction  convention  is  based  on  the 
SAE  (Society  of  Automotive  Engineers)  standard 
definitions  for  vehicle  dynamics.  Since  the  tire- 
terrain  models  are  currently  formulated  for  straight¬ 
ahead  rolling  (zero  slip  angle),  the  tire  longitudinal 
axis  is  in  the  direction  of  travel  and  the  lateral  or 
transverse  direction  is  perpendicular  to  travel  (in  the 
horizontal  plane),  as  shown  in  Eigure  6. 

The  tires  modeled  are  modem  radial  tires,  which  are 
constmcted  with  numerous  components,  as  shown  in  Eig¬ 
ure  7.  Each  of  these  elements  contributes  to  the  stmctural 
behavior  of  the  tire  and  is  considered  either  individually 
or  en  masse  when  creating  deformable  tire  models. 


Figure  6.  Tire  direction  convention. 
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Figure  7.  Construction  components  of  a  radiai  piy  tire. 
(After  Goodseii  1995.) 


Tire  deflection  is  a  primary  measure  of  the  tire 
structural  response  to  load  and  is  often  reported  based 
on  inflation  pressure  and  vertical  load.  Deflection  is 
the  difference  between  the  unloaded  and  the  loaded 
section  height  (and  is  sometimes  reported  in  percent¬ 
age,  normalized  by  the  unloaded  section  height).  The 
tire  cross  section  and  the  deformation  of  the  cross  sec¬ 
tion  under  load  are  illustrated  in  Figure  8.  A  thorough 
discussion  of  tire  mechanics  is  given  in  Clark  (1981). 


Figure  8.  Measurement  of  tire  dimensions. 


TERRAIN  MATERIAL  MODELS 

Terrain  substrate  subjected  to  vehicle  loads  has 
been  represented  using  a  wide  variety  of  material 
models,  including  elastic,  non-linear  elastic,  visco¬ 
elastic  (Pi  1988),  and  elastic-viscoplastic  (Saliba 
1990).  Recent  studies  have  concentrated  on  using 
either  cap  Drucker-Prager  plasticity  (Aubel  1993, 
1994,  Fervers  1994)  or  critical  state  models  such  as 
the  Bailey  and  Johnson  (1989)  soil  compaction  model 
implemented  by  Foster  et  al.  (1995)  or  a  new  critical 
state  model  (similar  to  Lade  and  Kim  1995)  imple¬ 
mented  by  Liu  and  Wong  (1996).  Comparisons  be¬ 
tween  the  Drucker-Prager  and  Cam-clay  models  for 
tire-terrain  interaction  have  been  published  by 
Meschke  et  al.  (1996)  for  snow  and  by  Chi  and 
Tessier  (1995)  for  soil.  Although  the  Cam-clay 
model  is  perhaps  the  most  widely  used  soil  model, 
both  of  the  above  studies  comment  on  convergence 
problems  when  using  Cam-clay.  The  choice  of  mate¬ 
rial  model  is  based  on  balancing  the  type  of  behavior 
desired  in  the  model  with  the  information  available  to 
determine  model  parameters.  The  terrain  materials 
simulated  in  this  study  were  fresh  snow  and  com¬ 
pacted  sandy  soil.  The  two  material  models  used 
were  a  modified  Drucker-Prager  with  a  cap,  and  a 
critical  state,  crushable  foam  model.  Both  models 
were  considered  suitable  for  capturing  the  highly 
compressible  behavior  of  natural  fresh  snow.  The 
compacted  sand  was  modeled  using  the  modified 
Drucker-Prager  cap  only.  A  description  of  these 
models,  model  parameter  determination,  and  valida¬ 
tion  of  the  material  models  with  test  data  follows. 

General  concepts  for  plasticity  models 

The  purpose  of  a  plasticity  model  is  to  define  the 
permanent  deformation  of  a  plastic  material.  This  is 
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done  with  the  following  basic  components  (Wood  and  Goodier  1970;  Jaeger  and  Cook  1969),  the  stress 
1990):  invariants  are: 


•  Elastic  properties  to  define  the  recoverable  de¬ 
formation; 

•  A  mathematical  surface  to  define  the  yield 
boundary  between  elastic  behavior  and  plastic 
material  behavior; 

•  A  plastic  flow  potential  to  mathematically  de¬ 
fine  the  plastic  deformation  (also  called  a  plas¬ 
tic  flow  law);  and 

•  A  hardening/softening  rule  defining  the  move¬ 
ment  (expansion  or  contraction)  of  the  yield 
surface  during  plastic  deformation. 

A  good  overview  of  the  development  of  plasticity 
theory  and  constitutive  modeling  of  soil  is  given  in 
Scott  (1985).  Schofield  (Schofield  and  Wroth  1968) 
extended  plasticity  theory  to  the  critical  state  concept, 
defining  either  contractile  or  dilatant  deformation  of 
porous  material  as  a  function  of  its  specific  volume  or 
void  ratio.  In  critical  state  theory  this  rule  is  devel¬ 
oped  around  the  concept  of  a  “critical  state,”  where 
the  plastic  shearing  deformation  occurs  at  a  constant 
volume.  Perhaps  the  most  famous  critical  state 
model,  the  Cam-clay  model,  was  developed  based  on 
the  behavior  of  clays.  The  concepts  are  equally  appli¬ 
cable  to  defining  the  shearing  and  volumetric  behav¬ 
ior  of  granular  materials  such  as  granular  soils  or 
snow  (Wood  1990).  Although  the  concepts  are  appli¬ 
cable  for  both  cohesive  and  granular  materials,  the 
behavior  of  the  granular  materials  has  not  been  ex¬ 
plored  as  thoroughly,  particularly  regarding  the  influ¬ 
ence  of  deviatoric  stress  on  the  yield  surface,  which 
is  less  clearly  defined  in  soils  but  may  take  on  a  much 
different  shape  than  the  yield  surface  of  metals. 
Wood  (1990)  considers  this  a  difference  in  detail 
rather  than  a  difference  in  concept. 

The  critical  state  class  of  model  used  in  this  study  is 
the  crushable  foam  model,  specifically  designed  for 
highly  compressible  materials,  is  which  a  characteristic 
of  fresh  snow.  The  modified  Drucker-Prager  cap 
model  also  has  the  features  of  a  critical  state  model 
(i.e.  regions  of  constant  volume  shear  deformation,  and 
compactive-dilatant  flow).  Both  models  use  non- 
associative  flow  (i.e.  the  flow  potential  is  not  associ¬ 
ated  with  the  yield  surface),  except  on  the  cap  surface 
of  the  Drucker-Prager  model. 

Yield  surface 

The  yield  surface  for  both  material  models  is  de¬ 
scribed  in  terms  of  stress  invariant  functions  of  the 
normal  Oi  and  shear  stress  Xy,  where  i  and  j  are  direc¬ 
tions  X,  y,  and  z  and  represent  principal  stresses  when 
i  =  1,  2,  or  3.  Based  on  elasticity  theory  (Timoshenko 


li  =  invariants  of  normal  stress: 
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where  the  components  of  stress  deviation  are 

=  Gj  -  p  for  i  =  X,  y,  z,  and  i  =  1,  2,  3.  (7) 

Model  parameters  are  defined  in  either  the  pres- 
sure-deviatoric  plane  (also  called  the  meridinal,  p-q, 
or  p-t  plane)  or  the  pressure-volume  plane  (v-p  or  v- 
In  p  plane)  according  to  the  following  definitions. 

1.  Mean  (total)  normal  stress,  also  called  equiva¬ 
lent  pressure  stress,  or  the  octahedral  normal  stress, 
which  determines  uniform  compression  or  dilation: 

P=<^oa=lh-  (8) 
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For  soils,  p  is  often  taken  as  the  mean  effective  stress 
p\  taking  into  account  the  pore  water  pressure,  u\ 

p'=p-u.  (9) 

2.  Deviatoric  stress,  also  called  the  Mises  equiva¬ 
lent  stress,  or  the  octahedral  shear  stress,  which  de¬ 
termines  distortion: 

^  =  %c^=^l|(h^+^h)=S^  (10) 

3.  An  additional  deviatoric  stress  measure,  t  (HKS 
1998): 


,  1 

(  T  ^ 

3" 

t  =  — 

1+ — 

2 

K 

1  kJ 

1  ^  J 

where  ^  is  a  material  parameter  (between  0.778  and 
1.0)  that  controls  the  yield  dependence  on  the  third 
stress  invariant  and  defines  the  shape  of  the  yield 
surface  in  the  deviatoric  plane.  The  K  value  allows 
for  different  stress  failures  in  tension  and  compres¬ 
sion.  For  K  =  I,  the  yield  surface  is  circular  (von 
Mises  yield),  as  shown  in  Figure  9,  and  the  failure 
stress  is  the  same  in  tension  and  compression.  K  is 
limited  to  values  of  0.778  or  greater  for  the  surface  to 
remain  convex.  Some  of  the  more  common  failure 
surfaces  are  shown  in  Figure  10  in  deviatoric  space 
and  in  Figure  1 1  in  principal  stress  space. 


Figure  9.  Modified  Drucker-Prager  yieid  sur¬ 
face  in  deviatoric  space.  (After  HKS  1998.) 


S3 


Figure  10.  Yieid  surfaces  in  deviatoric  space. 
(After  HKS  1998.) 


Figure  11.  Common  yieid  surfaces  aiong  the  hy¬ 
drostatic  axis  in  principai  stress  space.  Caps 
defining  compactive  faiiure  are  not  shown.  (After 
Shen  and  Kushwaha  1998.) 


7 


to  contents 


Cap  plasticity  models  incorporate  both  the  com¬ 
pression  and  shearing  of  the  material.  Critical  state 
models  are  a  type  of  cap  model  where  a  critical  state 
line  (csl)  is  located  at  the  peak  of  the  yield  surface 
and  divides  material  yield  by  compaction  on  one  side 
and  yield  by  dilation  on  the  opposite  side.  The  mate¬ 
rial  deforms  at  a  constant  volume  at  the  critical  state. 
The  critical  state  line  (curve)  in  three-dimensional 
space  defined  by  normal  compressive  stress  p  \  devia- 
toric  (shear)  stress  q,  and  specific  volume  v  is  shown 
in  Figure  12.  For  simplicity  in  gathering  the  model 
parameters,  the  yield  surface  is  generally  viewed  in 
planar  form  in  the  compressive  plane  (v-p  or  v-ln  p) 
or  the  deviatoric  plane  (p'-q  or  p-t). 

Hardening  law 

The  cap  on  the  yield  surface  defines  the  compac- 
tive-dilatant  behavior  of  the  material  failure  as  the 
cap  expands  and  contracts  (Fig.  13).  The  cap  is  gen¬ 
erally  spherical  or  ellipsoidal,  and  the  material  either 
hardens  or  softens  by  expanding  or  contracting  the 
cap.  This  behavior  is  defined  in  the  pressure- volume 
relationship  called  a  hardening  law.  The  pressure- 
volume  relationship  is  often  exponential  and  there¬ 
fore  can  be  modeled  using  an  exponential  hardening 


law  (Fig.  14).  In  cases  where  exponential  hardening 
is  not  a  good  fit,  as  is  common  in  soils,  the  hardening 
law  can  also  be  represented  in  a  piecewise  linear  ap¬ 
proach  using  the  experimental  data  (Fig.  15). 


Figure  12.  Three-dimensional  view  of  critical  state 
yield  surface.  The  critical  state  line  defines  failure 
at  constant  volume.  (After  Wood  1990.) 


Figure  13.  Cap  contraction  or  expansion  reflecting  the  softening  (left) 
or  hardening  (right)  of  the  material.  (After  HKS  1998.) 
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Figure  14.  Exponential  hardening  law  for  material  in  hydro¬ 
static  compression. 


Figure  15.  Piecewise  linear  modeling  of  the 
hardening  law.  (After  HKS  1998.) 


Hardening  can  be  modeled  as  isotropic  hardening, 
where  the  yield  surface  grows  equally  in  all  direc¬ 
tions,  or  as  kinematic  hardening,  where  the  yield  sur¬ 
face  moves  as  a  rigid  body  in  principal  stress  space 
(Fig.  16).  Kinematic  hardening  is  more  accurate  for 
modeling  cyclic  loading  of  soils  but  is  more  compu¬ 
tationally  intensive.  A  combination  of  the  two  types 
of  hardening  is  probably  the  best  representation  for 
cyclic  loading  of  soil.  The  Multi-Mechanical  Model 
(MMM)  proposed  by  Peters  (1991)*  implements 


*  Also  personal  communication.  Geotechnical  and 
Structures  Laboratory,  ERDC,  Vicksburg,  MS,  1999. 


some  of  these  features  without  the  cumbersome 
tracking  of  yield  surfaces  by  engaging  a  series  of 
yield  surfaces  (failure  mechanisms).  The  MMM  con¬ 
cept  has  been  successfully  applied  to  modeling  the 
cyclic  loading  of  pavements  by  Smith  (2000).  This 
constitutive  model  is  still  under  development  and  is 
currently  being  calibrated  for  loose,  thawing  soil  be¬ 
havior  for  future  applications  to  tire-terrain  model¬ 
ing.  The  material  models  implemented  in  this  study 
use  isotropic  hardening,  which  is  considered  a  rea¬ 
sonable  approximation  for  monotonic  loading  (one 
pass  of  the  wheel). 
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M  =  slope  of  the  critical  state  line  in  the  p-q 
plane,  which  controls  the  shape  of  the 
yield  ellipse 

E  =  Young’s  modulus  of  elasticity  defining 
recoverable  strain 
V  =  Poisson’s  ratio. 

The  following  sections  explain  the  Drucker- 
Prager  cap  model  and  the  crushable  foam  model  in 
more  detail. 

Modified  Drucker-Prager  cap  model 

For  the  modified  Drucker-Prager  cap  model  used 
in  this  study,  the  yield  surface  is  a  modified  von 
Mises  yield  (i.e.,  the  material  constant  K  =  1.0), 
which  is  circular  in  the  deviatoric  plane  (Fig.  9).  In 
the  p-t  plane  the  yield  surface  has  two  major  seg¬ 
ments:  the  Drucker-Prager  shear  portion  of  the  curve 
(analogous  to  the  Mohr-Coulomb  line)  defines  shear 
deformation,  and  the  cap  portion  of  the  surface  de¬ 
fines  the  intersection  with  the  pressure  axis  (Fig.  17). 
The  following  equations  define  the  yield  criteria  in 
each  section  of  the  yield  surface. 

For  the  Drucker-Prager  shear  or  distortional  fail¬ 
ure: 


Fg  =  r - /? tanP- J  =  0  .  (12) 

Kinematic 


Figure  16.  Isotropic  and  kine¬ 
matic  hardening.  (After  Shen 
and  Kushwaha  1998.) 

Parameters 

A  relatively  small  number  of  parameters  are 
needed  for  these  material  models.  Although  more 
detailed  models  may  provide  a  more  accurate  solu¬ 
tion,  the  parameters  needed  to  define  the  model 
become  prohibitively  numerous  and  are  often  un¬ 
available.  The  objective  of  any  model  is  to  capture 
the  most  important  material  behaviors  without  being 
unduly  cumbersome.  The  following  descriptions  of 
the  material  parameters  are  based  on  the  models  as 
implemented  in  ABAQUS  (HKS  1998).  Additional 
information  on  soil  plasticity  and  critical  state  soil 
mechanics  is  given  in  Wood  (1990). 

The  model  parameters  used  in  this  study  are: 

d  =  Drucker-Prager  material  cohesion 
P  =  Drucker-Prager  material  angle  of  friction 
X  and  K  =  slopes  of  the  loading  (compression)  and 
unloading  (elastic)  lines  in  the  compression 
plane  as  indicated  in  Figure  14.  This  hard¬ 
ening  or  softening  of  the  material  can  also 
be  defined  in  a  piecewise  manner  (Fig.  15). 


For  the  cap  region  of  compactive-dilatant  failure: 


^c=J(p-Pa)  + 


Rt 


12 


(1-ha-a/cosP) 
-R(d  +  p^  tanP)  =  0 


(13) 


where  a  is  a  transition  parameter,  ranging  typically 
from  0.0  to  0.05,  that  smooths  the  transition  between 
the  shear  failure  and  the  cap  failure,  is  a  material 
parameter  controlling  the  cap  eccentricity,  and  Pa 
defines  the  evolution  of  the  cap  along  the  pressure 
axes  (via  /^b-^voi  hardening  law)  according  to 


P\y-Rd 
(l  +  7?tanp) 


(14) 


The  transition  between  the  shear  and  the  cap  fail¬ 
ure  is 


r=J(p-pJ  + 


t- 


1-- 


a 


cosP 


{d  +  p^  tanP) 


-a  (J  +  tan  p)  =  0  . 


(15) 
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Figure  17.  Modified  Drucker-Prager  cap  yieid  surface  in  the  p-t 
piane.  (After  HKS  1998.) 


The  pressure-volume  relationships  define  both 
hardening  and  softening  through  volume  changes 
based  on  how  the  cap  portion  of  the  yield  surfaces 
expands  and  contracts.  The  hardening  of  the  cap  can 
be  defined  as  an  exponential  relationship  in  the  pres¬ 
sure-volume  space  as  shown  in  Figure  14,  or  it  is 
entered  in  a  piecewise  fashion  (Fig.  15)  as  a  table  of 
and  eSoi  pairs.  The  piecewise  approach  is  recom¬ 
mended  by  HKS  (1998)  for  a  better  fit  to  the  data  and 
better  model  performance.  Both  methods  were  used 
in  this  study. 

The  plastic  flow  is  defined  by  an  elliptical 
potential  surface  (Fig.  18).  Flow  is  associative  (nor¬ 


mal  to  the  surface)  in  the  cap  region;  therefore,  the 
equation  for  the  flow  surface  is  identical  to  the  equa¬ 
tion  for  the  cap  yield  surface.  In  the  transition  and 
shear  regions,  the  flow  is  non-associative  (the  flow 
potential  is  not  perpendicular  to  the  failure  surface) 
and  is  defined  by  the  following  equation  for  the  flow 
potential: 


Gs  =  J[(Pa-p)tanpf  + 


n2 


1-ha-a/cosP 


(16) 


Figure  18.  Modified  Drucker-Prager  cap  fiow  potentiai  in  the  p-t 
piane.  (After  HKS  1998.) 


11 


to  contents 


Figure  19.  Yield  surface  for  the  crushable  foam  model. 
(After  HKS  1998.) 


Figure  20.  Plastic  flow  surface  for  the  crushable  foam  model. 


Crushable  foam  model 

The  crushable  foam  model  is  a  class  of  critical 
state  model  and  has  many  similarities  to  the  form  of 
the  Drucker-Prager  cap  model.  The  crushable  foam 
model,  however,  has  a  single  ellipse  forming  the 
yield  surface  (Fig.  19),  defined  by  (HKS  1998) 

^foam  ~  f  ~  fo  (17) 


/  .  \21 

LI  2  J 

-h 

1 

where  is  the  material  strength  in  hydrostatic  ten¬ 
sion,  Pc  is  the  material  yield  stress  in  hydrostatic 
compression,  and  M  is  the  slope  of  the  critical  state 
line  in  p-t  plane.  As  indicated  in  Figure  19,  the  yield 
surface  hardens  and  softens  with  yield  but  the  pt  is 
assumed  constant  (isotropic  hardening). 

In  the  crushable  foam  model  the  plastic  flow  is 
non-associative  and  perpendicular  to  the  flow  surface 
^  (Fig.  20): 

S  = 


This  assumes  a  flow  direction  identical  to  the 
stress  direction  for  stress  paths  along  the  p  axis, 
which  suggests  that  loading  in  any  principal  direction 
results  in  insignificant  deformation  in  the  other  direc¬ 
tions  (HKS  1998).  The  validity  of  this  assumption  for 
snow  was  examined  as  part  of  this  study. 

The  hardening  law  takes  the  same  form  as  the 
hardening  law  for  the  Drucker-Prager  cap  model. 

Constitutive  models  for  snow 

Two  material  models  were  judged  to  have  per¬ 
formance  characteristics  most  similar  to  natural  fresh 
snow:  a  Drucker-Prager  cap  plasticity  model  and  a 
critical  state  plasticity,  crushable  foam  model.  Liu 
(1994)  used  the  Drucker-Prager  plasticity  model  for 
modeling  the  sliding  of  rubber  blocks  on  snow. 
Mundl  et  al.  (1997)  extended  Liu’s  work  using  a 
multi-surface  plasticity  model  to  optimize  the  snow 
behavior  under  both  shear  and  compression.  Their 
intent,  however,  was  to  simulate  shearing  forces  on 
compacted  snow  on  roads  (with  a  density  of  500 
kg/m^),  which  is  a  significantly  different  material 
than  encountered  during  cross-country  mobility  on 
fresh  snow  (with  a  density  of  200-300  kg/m^).  For 
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this  study  the  plasticity  model  used  was  the  well- 
documented  Drucker-Prager  cap. 

The  second  material  model,  a  crushable  foam 
model,  was  chosen  based  on  work  by  Johnson  et  al. 
(1992),  where  a  crushable  foam  model  was  success¬ 
fully  applied  to  natural  snow.  Johnson  used 
PRONTO,  an  explicit  finite  element  code  from  San- 
dia  National  Lab,  to  simulate  shock  wave  propaga¬ 
tion  in  one  dimension.  For  tire-snow  interaction, 
however,  the  model  must  simulate  snow  deformation 
in  a  three-dimensional  stress  field.  Initial  simulations 
using  the  ABAQUS  crushable  foam  model  showed 
that  no  lateral  deformation  resulted  from  an  axial 
load.  This  issue  prompted  questions  regarding  the 
deformation  behavior  of  snow  under  a  purely  vertical 
load  and  whether  the  predicted  response  from  a 
crushable  foam  model  was  suitable  for  snow  in  a 
three-dimensional  stress  state.  Experiments  measur¬ 
ing  the  lateral  deformation  of  snow  under  a  purely 
vertical  load  were  conducted  in  cooperation  with  the 
Keweenaw  Research  Center  in  Houghton,  MI  (Shoop 
and  Alger  1998).  Results  indicate  that  vertical  load¬ 
ing  of  non-sintered  snow  causes  primarily  vertical 
deformation.  Fukue  (1979)  performed  experiments 
that  showed  that  lateral  deformation  was  minimal  for 
high  load  rates.  Both  of  these  studies  support  our 
choice  of  crushable  foam  models  for  vehicle  loading. 

Material  parameter  determination 

Extreme  changes  in  snow  properties  with  time,  tem¬ 
perature  gradients,  applied  load,  and  deformation  pro¬ 
hibit  model  parameter  acquisition  from  the  same  snow. 
Therefore,  approximations  were  made  by  estimating 
parameters  using  test  data  from  snow  of  similar  charac¬ 
teristics  (density,  age,  snow  type,  and  location).  The 
type  of  snow  modeled  was  a  fresh  snow  with  a  density 
of  200-250  kg/m^  at  moderate  temperatures  (between 
-10°  and  -1°C).  Test  data  were  gathered  from  the  field 
and  from  the  literature  to  match  this  snow  type  as 
closely  as  possible.  Because  data  were  gathered  from 
several  studies,  the  snow  was  not  exactly  the  same  in 
crystal  structure,  but  this  effect  is  minimized  in  fresh 
snow  where  sintering  (bonding  between  snow  grains) 
has  not  occurred.  Also,  most  of  the  tests  were  performed 
on  “lake  effect”  snow  from  the  Keewenaw  Peninsula  in 
the  Upper  Peninsula  of  Michigan,  which  is  a  fairly  con¬ 
sistent  snow  type  under  normal  winter  conditions.  A 
discussion  of  selecting  the  initial  values  of  the  material 
parameters  follows. 

Elastic  properties.  For  snow  deformation  under  a 
tire,  the  plastic  deformation  is  much  greater  than  the 
elastic  deformation;  however,  the  elastic  deformation 
contributes  to  energy  losses  during  rolling  resistance. 
The  basic  elastic  parameters.  Young’s  modulus  E  and 


Poisson’s  ratio  v,  were  estimated  based  on  data  com¬ 
piled  by  Mellor  (1975)  and  expanded  on  by  Shapiro 
et  al.  (1997)  (Fig.  21).  Since  the  elastic  modulus  is  a 
strong  function  of  density,  varying  over  four  orders 
of  magnitude,  the  elastic  modulus  should  ideally  be 
modeled  as  a  function  of  the  density  as  a  state  vari¬ 
able.  For  this  study,  however.  Young’s  modulus  was 
held  constant  at  a  value  equivalent  to  snow  with  a 
density  of  300-400  kg/m^.  At  snow  densities  less 
than  this,  the  elastic  contribution  will  be  minimal. 


A  =  Pulse  propagation  or  flexural  vibration  at  high 
frequencies,  -10°  to  -25  °C 

B  =  Uniaxial  compression,  strain  rate  approximately  3 
X  10“^  to  2  X  10“^  s“\  temperature  -25°C 
Ci  =  Uniaxial  compression  and  tension,  strain  rate  ap¬ 
proximately  8  X  10“^  to  4  X  10^  s“\  -12°  to  -15°C 
C2  =  Static  creep  test,  -6.5°  to  -19°C 
D  =  Complex  modulus,  10^  Hz,  -14°C 
K  =  Static  Young’s  modulus  and  quasi-static  Pois¬ 
son’s  ratio 

S  =  Quasi-static  measurements  of  Poisson’s  ratio 

Figure  21.  Compilation  of  Young’s  modulus  and 
Poisson’s  ratio  measurements  on  snow.  [After 
Mellor  (1975)  and  Shapiro  et  al.  (1997)] 
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Porous  elasticity  accounts  for  the  elastic  properties 
of  materials  that  undergo  large  volume  changes,  such 
as  fresh  snow.  To  implement  a  porous  elasticity 
model,  the  additional  parameters  needed  are  the  slope 
of  the  elastic  portion  of  the  pressure-volume  curve 
(during  an  unload-reload  event,  as  illustrated  in  Fig¬ 
ure  14);  K,  representing  the  elastic  portion  of  compac¬ 
tion;  and  the  elastic  tensile  limit,  pf .  The  elastic 
compaction  value  was  not  available  since  compaction 
tests  were  performed  under  monotonic  loading  (no 
unloading),  so  k  was  chosen  to  provide  realistic  elas¬ 
tic  behavior  within  the  stress  ranges  of  the  simula¬ 
tions.  The  elastic  tensile  limit  pf  was  chosen  to  be 
the  same  as  the  same  as  tensile  limit  p^  used  to  define 
the  yield  surface. 

Yield  surface.  The  parameters  defining  the  shear 
portion  of  the  yield  surface  for  the  Drucker-Prager 
model  can  be  calculated  from  the  Mohr-Coulomb 
cohesion  c  and  the  internal  angle  of  friction  (|),  assum¬ 
ing  plane  strain  response  and  non-dilatant  flow.  Al¬ 
though  non-dilatant  flow  can  be  assumed  for  the 


compaction  beneath  the  tire,  a  plane  strain  response 
may  be  more  appropriate  for  deep  snow  than  for  shal¬ 
low  snow,  as  illustrated  in  the  snow  deformation 
measurements  in  Figure  3.  The  Drucker-Prager  mate¬ 
rial  cohesion  d  and  the  Drucker-Prager  material  an¬ 
gle  of  friction  (3  can  be  calculated  from 

tan  (3  =  1.73  sin  (|)  (19) 

d  =  c  1.73  cos  (|).  (20) 

Acquiring  apparent  values  for  cohesion  c'  and  fric¬ 
tion  angle  (|)'  from  shear  loading  of  snow  under  a  ve¬ 
hicle  tire  and  from  a  ring  shear  device  is  discussed  by 
Blaisdell  et  al.  (1990)  and  Alger  and  Osborne  (1989). 
Both  methods  were  applied  to  field  measurements  on 
undisturbed  snow  with  a  density  ranging  from  60  to 
250  kg/no?  and  temperatures  from  -2°  to  -16°C.  The 
results  are  shown  in  Figure  22,  along  with  the  values 
used  in  the  CRREL  Shallow  Snow  Model  (SSM  1.0) 
(Blaisdell  et  al.  1990). 


Figure  22.  Calculated  c'  and  using  measurements 
from  the  CRREL  Instrumented  Vehicle  (CIV)  and  the 
ring  shear  device  and  predicted  using  the  Shallow 
Snow  Model,  SSM  1.0.  (From  Blaisdell  et  al.  1990.) 
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Using  these  graphs,  for  a  snow  density  of  200 
kg/m^,  values  ranged  from  c'  =  2.14  kPa  and  (|)'  = 
8.98°  for  the  ring  shear  device  to  c'  =  43  kPa  and  (|)'  = 
14°  using  the  CIV.  Because  of  the  nature  of  the  test, 
these  values  may  be  more  representative  of  the  inter¬ 
face  shear.  However,  similar  values  of  cohesion  (ap¬ 
proximately  20-30  kPa)  were  reported  in  Shapiro  et 
al.  (1997).  These  values  yield  a  Drucker-Prager  (3 
ranging  from  15.2°  to  22.5°  and  a  d  ranging  from  3.7 
to  72  kPa. 

Parameters  describing  the  shape  of  the  cap  were 
adjusted  based  on  model  response.  The  flow  stress 
ratio  K  (ratio  of  tensile  strength  to  compressive 
strength)  was  held  at  1.0.  This  agrees  with  data  pre¬ 
sented  in  Mellor  (1975)  and  Shapiro  et  al.  (1997), 
indicating  nearly  equal  values  of  compressive  and 
tensile  strength  for  low-density  snow  (data  set  M  in 


Fig.  23).  The  cap  eccentricity  parameter  R  was  cho¬ 
sen  based  on  typical  values  for  earth  materials  having 
a  very  steep  compaction  cap  {R  =  0.02  to  0.0001). 
The  transition  surface  radius  a  was  set  to  0.0  for  no 
transition.  The  initial  cap  yield  surface  position  eJoJo 
was  arbitrarily  set  to  0.001  to  allow  initial  softening. 

For  the  foam  model  the  initial  yield  stress  in  uni¬ 
axial  compression  Gq  was  based  on  data  in  Shapiro  et 
al.  (1997).  The  initial  yield  surface  position  in  hydro¬ 
static  compression  Pcio  (for  exponential  hardening) 
was  chosen  to  allow  for  initial  softening  or  hardening 
(about  half  the  value  of  gq).  The  hydrostatic  tensile 
strength  was  selected  as  10%  of  /?cio  as  a  general 
rule  for  foam  materials.  M  is  calculated  from  the 
other  input  parameters. 

These  material  parameters  were  refined  based  on 
model  response  and  are  presented  in  Table  1. 


Table  1.  Initial  material  model  parameters. 

Drucker-Prager  cap 

Foam 

Elastic 

of  =5, 000  to  30,000  Pa 

A  =  0.35 

K=  0.005 

15  =  22.53° 

Pcio=  113.76  Pa 

v  =  0.3 

R=2.2E-2  to  1.1  E-4 

-e?oilo+e?|„r  =0.0005 

pf  =  11 .376  Pa 

=0.001 

Pclo+  pt=  1 13.76  Pa 

E=  13.79  MPa 

a=  0.0 

Pt=  11.376  Pa 

K=  1.0 

(70=275.79  Pa 

K=  1.0 

A  =  Strain  rates  greater  than  5x10“"^  s“^ 

B  =  In- situ  tensile  strength 
M  =  Uniaxial  compressive  and  tensile  strength 
under  rapid  loading 

Figure  23.  Compilation  of  uniaxial  strength  data. 
[After  Mellor  (1975)  and  Shapiro  et  al.  (1997).] 


Hardening.  Two  methods  were  used  to  handle  the 
hardening  data:  1)  using  an  exponential  hardening 
function  and  2)  using  tabular  data.  For  exponential 
hardening,  the  slope  of  pure  hydrostatic  compression 
behavior  X  is  obtained  by  plotting  the  void  ratio 
^  against  ln[(p  Pi)/pi],  as  shown  in  Figure  14,  where  p 

i  is  the  hydrostatic  pressure  and  p^  is  the  strength  in 

«  hydrostatic  tension.  For  the  snow  pressure-volume 
w  relationship,  void  ratio  was  plotted  against  ln(p)  since 
p^  was  not  precisely  known.  For  snow,  p^  is  a  very 
small  quantity  and  only  the  slope  is  required,  which 
is  minimally  affected  when  p»Pi.  These  plots  were 
generated  using  the  compression  test  data  published 
in  Abele  and  Gow  (1975). 

Values  for  the  logarithmic  hardening  modulus 
were  calculated  using  three  different  compression 
tests  of  similar  snow  density,  along  with  values  cal¬ 
culated  from  the  final  pressure-density  pairs  of  an 
additional  21  tests,  all  from  Abele  and  Gow  (1975). 
Since  the  measurements  are  from  one-dimensional 
compression  or  consolidation  (oedometer)  tests  by 
application  of  a  vertical  load  to  a  sample  within  a 
rigid  cylinder  (i.e.  uniaxial  strain),  the  calculation  of 
the  mean  hydrostatic  pressure  p  is  not  straightfor¬ 
ward.  One  approximation  calculates  the  horizontal 
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(24) 


stresses  Gh  based  on  earth  pressure  theory  assuming 
that  they  are  related  to  the  vertical  stress  Gy  by  a  fac¬ 
tor  of  ^0^  known  in  soil  mechanics  as  the  earth  pres¬ 
sure  coefficient,  where 

C!t=KoO^.  (21) 

The  earth  pressure  factor  can  be  calculated  from 
the  friction  angle  using 

=  1  -  sin  (p .  (22) 

Based  on  Alger’s  shear  measurements  (Alger 
1988,  Alger  and  Osborne  1989),  the  average  value 
is  approximately  0.12.  Assuming  that  the  maximum 
principal  stress  is  vertical  (as  would  be  expected  in 
uniaxial  compression  tests  on  undisturbed  snow)  and 
using  equations  8  and  21: 

p  =  \{Gi+a2+c^)  =  \{l  +  2Ko)Gi  (23) 

which,  for  =  0.12,  reduces  to 


p  =  0.413ai . 

Using  this  expression  to  calculate  the  mean  stress 
from  individual  uniaxial  consolidation  tests  with  a 
density  of  200-390  kgW  (Abele  and  Gow  1975, 
tests  3,  10,  and  29)  results  in  a  range  of  X  between 
0.223  and  0.454. 

An  additional  data  set  was  obtained  using  the 
combined  data  from  21  consolidation  (oedometer) 
tests  on  snow  with  initial  densities  of  200-220  kg/no? 
and  temperatures  of  0°  to  -3°C.  These  data  are  shown 
in  Figure  24  and  fall  within  the  expected  range  of 
stress-density  data  for  snow  shown  in  Figure  25. 
Selected  points  from  the  best-fit  line  of  these  data 
were  used  to  calculate  the  values  in  Table  2  and  are 
plotted  in  Figure  26  to  obtain  an  average  A.  of  0.33. 

The  hardening  table  data  were  taken  from  an  aver¬ 
age  of  Abele  and  Gow  tests  on  warm,  200-kg/m^ 
snow  (tests  10  and  29)  and  are  given  in  Table  3. 


Figure  24.  Compression  of  naturai  snow  at  0°  to  -3°C.  (After  Abeie 
and  Gow  1975.) 


Tabie  2.  Pressure-voiume  data 

caicuiated  from  the 

Abeie  and 

Gow  (1975)  data  set  for  200- 

to  220-kg/m^ 

snow  at  0° 

to  -3X. 

P 

Oi 

p  =  0.413  Oi 

Void  ratio, 

(kg/m") 

(MPa) 

(MPa) 

In  p 

e 

900 

2.8 

1.16 

0.15 

0.00 

770 

2 

0.83 

0.19 

0.17 

620 

1 

0.41 

0.88 

0.45 

540 

0.5 

0.21 

-1.58 

0.67 

470 

0.2 

0.08 

-2.49 

0.91 

410 

0.1 

0.04 

-3.19 

1.20 

380 

0.05 

0.02 

-3.88 

1.37 
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a^,  Major  Principal  Stress  (MPa) 


Figure  25.  Compression  of  undisturbed 


A  =  Natural  densification  of  snow  at  -1°  to  -48°C 
B  =  Slow  natural  compression  of  dense  firn  and 
porous  ice  (from  polar  caps) 

C  =  Slow  compression  of  solid  ice 
E  =  Calculated  values  from  plane  wave  impact  at 
20-40  m/s 

F  =  Hugoniot  data  for  explosively  generated  shock 
waves  with  impact  velocity  of  1-12  m/s  at  tem¬ 
peratures  of  -7°  to  -18°C 

J  =  Compression  at  strain  rate  10“^  s“^  at  -7°  to 
-18°C 

K  =  Compression  in  uniaxial  strain  with  incre¬ 
mental  loading  at  -2°  to  -3°C 
Solid  Lines  =  Abele  and  Gow  data,  T  =  0°  to  -34°C, 
po  =  90  to  270  kg/m^ 

w  and  ice.  (After  Abeie  and  Gow  1975.) 


Figure  26.  Pressure-voiume  data  from  compression  tests  of  Abeie  and  Gow  (1975). 
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Table  3.  Hardening  table  for  fresh  snow  with  a  density  of 
200  kg/m^. 


p  =  1/3  Oi  (Pa) 

p  =  0.413  Oi  (Pa) 

p  =  ai  (Pa) 

pPi 

^vol 

37.92 

46.98 

113.76 

0 

0.017E6 

0.021  E6 

0.05E6 

0.593 

0.033E6 

0.041  E6 

0.1  E6 

0.669 

0.067E6 

0.083E6 

0.2E6 

0.806 

0.167E6 

0.207E6 

0.5E6 

0.944 

0.333E6 

0.413E6 

1.0E6 

1.083 

0.667E6 

0.826E6 

2.0E6 

1.299 

0.933E6 

1.156E6 

2.8E6 

1.455 

1.083E6 

1 .342E6 

3.25E6 

1.475 

2.0E6 

2.478E6 

6.0E6 

1.50 

2.0E7 

2.478E7 

6.0E7 

1.514 

Figure  27.  Comparison  of  model  and  experimental  data  for  uniaxial 
compression  tests  on  low-density  snow. 


Model  calibration  against  laboratory  strength 
tests.  A  finite  element  model  simulating  the  uniaxial 
confined  compression  tests  was  used  to  refine  the 
material  properties  shown  in  Tables  1  and  3.  A  single 
element  model  was  determined  to  be  adequate  after 
multi-element  models  gave  the  same  results.  The 
uniaxial  compressive  test  data  used  for  calibration 
were  from  the  best-fit  line  for  fresh  snow  of  density 
200-220  kg/m^  and  temperatures  of  0°  to  -3°C  from 
Figure  24.  A  comparison  between  the  experimental 
data  and  model  results  is  given  in  Figure  27.  The 
agreement  between  modeled  and  measured  material 
behavior  is  good  for  both  the  Drucker-Prager  cap 
model  and  the  crushable  foam  model. 

Material  model  validation  with 
plate  sinkage  tests 

The  initial  model  parameters  were  applied  to  a 
plate  sinkage  test  for  snow  of  similar  age  and  density. 


Two  different  plate  sinkage  tests  were  used  in  the 
comparisons.  The  field  test  consisted  of  pushing  a  20- 
cm-  (8-in.-)  diameter  rigid  plate,  originally  at  the 
surface,  into  40  cm  (16  in.)  of  fresh  snow  with  aver¬ 
age  density  of  200  kg/m^  (Alger  and  Osborne  1989). 
The  laboratory  test  pushed  a  23 -cm-  (9-in.-)  diameter 
rigid  plate  into  a  cube  of  snow  60  x  60  x  50  cm  deep 
with  an  average  initial  density  just  under  200  kg/m^. 
Three  different  models  were  used  to  simulate  the 
plate  sinkage  tests:  1)  a  Drucker-Prager  cap  material, 
2)  a  Drucker-Prager  cap  material  with  vertical  shear 
surfaces  along  the  snow-plate  interface,  and  3)  a 
crushable  foam  material  model  with  shear  surfaces 
along  snow-plate  interface.  The  vertical  frictional 
surfaces  allowed  these  elements  to  slide,  reducing 
element  distortion,  which  more  closely  mimicked  the 
shear  surfaces  observed  in  the  field.  The  slide  planes 
were  modeled  using  a  contact  surface  with  Coulomb 
friction  of  jl  =  0.3. 
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b.  Drucker-Prager  cap  with  vertical  shear  contact 
surface. 


c.  Crushable  foam  with  shear  contact  surface. 

Figure  28.  Deformed  meshes  for  field  plate  sink- 
age  simulations. 


The  plate  was  lowered  by  constraining  the  surface 
nodes  radially  while  displacing  them  into  the  snow, 
effectively  creating  a  no- slip  contact  between  the 


plate  and  snow.  The  deformed  meshes  for  each  of  the 
models  of  the  field  test  are  shown  in  Figure  28.  Fig¬ 
ure  28b  clearly  shows  the  effect  of  the  shear  surface 
on  the  sides  of  the  plate,  with  the  model  deforming 
similarly  to  field  observations,  compared  to  the  unde¬ 
sirable  element  distortion  along  the  sides  of  the  plate 
when  these  elements  are  not  allowed  to  slide  (Figure 
28a).  In  future  models,  adaptive  meshing,  where  the 
mesh  adjusts  to  accommodate  large  distortions, 
should  help  alleviate  the  problem  of  excessive  distor¬ 
tion  without  the  use  of  the  shear  contact  surface.  Also 
clearly  seen  in  Figure  28  is  the  lateral  deformation  of 
the  snow  beneath  the  plate  in  the  Drucker-Prager  cap 
models  (Fig.  28a  and  b),  but  no  lateral  deformation  is 
seen  in  the  crushable  foam  model  (Fig.  28c). 

The  force-displacement  curves  for  the  plate  mov¬ 
ing  into  the  snow  are  compared  to  the  experimental 
measurements  in  the  laboratory  and  field  in  Figure 
29.  The  model  matches  the  controlled  laboratory  test 
well,  particularly  for  the  Drucker-Prager  cap  model. 
In  the  field  test  the  hardening  behavior  was  difficult 
to  simulate.  The  laboratory  test  was  not  run  to  maxi¬ 
mum  resistance  and  would  be  similar  to  a  “deep 
snow”  condition,  while  the  plate  in  the  field  test  was 
pushed  to  ultimate  resistance,  causing  compaction 
equivalent  to  a  “shallow”  snow  condition. 

Snow  density  can  be  calculated  from  the  modeled 
volumetric  plastic  strain  using  the  hardening  tables 
and  the  pressure-volume  relationships  in  Tables  2 
and  3.  The  snow  density  values  from  the  models  are 
shown  in  Figure  30.  These  densities  can  be  compared 
to  snow  densities  measured  under  the  plate,  as  shown 
in  Figure  31.  The  maximum  volumetric  inelastic 
strain  in  the  plate  sinkage  model  occurred  in  isolated 
locations  adjacent  to  the  plate  at  values  of  approxi¬ 
mately  0.5,  which  is  equivalent  to  a  snow  density  of 
420  kg/m^,  whereas  very  dense  snow  (440-520 
kg/m^)  was  measured  up  to  5  cm  below  the  plate 
(Fig.  31).  These  values  correspond  to  a  strain  of  0.5 
m/m  or  greater  (assuming  primarily  vertical  deforma¬ 
tion)  on  the  uniaxial  stress-strain  model  calibration 
curve  (Fig.  27).  The  divergence  of  the  models  and 
data  at  strains  greater  than  0.5  m/m  (Fig.  26)  would 
account  for  density  discrepancies  in  the  plate  sinkage 
measurements  and  models  in  high  strain  areas.  On  the 
whole,  however,  most  of  the  strain  occurs  at  levels  of 
0.4  or  less,  so  the  model  provides  a  good  fit  to  the 
overall  plate  sinkage  behavior. 

The  laboratory  tests  were  also  used  to  evaluate  the 
lateral  deformation  of  the  snow,  since  field  observa¬ 
tions  show  lateral  deformation  under  a  vehicle  load. 
On  closer  examination  of  the  literature,  however,  few 
examples  were  available  where  the  snow  had  been 
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a.  Laboratory  tests. 


Figure  29.  Comparison  of  modeied  and  measured  piate  sinkage  test  for 
iaboratory  and  fieid  tests. 


marked  in  a  way  that  would  enable  the  observation  of 
lateral  movement. 

In  addition,  the  nature  of  loading  from  vehicles  is 
not  strictly  vertical.  Therefore,  an  experiment  was 
devised  to  mark  several  snow  samples  both  vertically 
and  horizontally  and  subject  them  to  purely  vertical 
loading.  The  snow  in  Figure  31  was  one  such  ex¬ 
periment,  clearly  indicating  no  lateral  deformation, 
since  the  vertical  marks  remained  vertical  throughout 


the  deformation  (Shoop  and  Alger  1998).  Fukue 
(1979)  reported  similar  observations  in  light  snow 
(120  kg/m^)  as  a  function  of  penetration  speed,  as 
shown  in  Figure  32.  This  substantiates  the  use  of  the 
crushable  foam  model  to  represent  fresh  snow.  How¬ 
ever,  dense  or  bonded  snow  does  show  lateral  defor¬ 
mation  (Shoop  and  Alger  1998),  as  does  snow  sub¬ 
jected  to  slow  loading  (Fukue  1979)  (Fig.  32a). 
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a.  Laboratory  test  simulation. 


Figure  30.  Modeled  snow  density  (kg/m^)  for  plate  sinkage  test  for  Drucker-Prager  cap  material. 


22.86  cm 


Figure  31.  Measured  displacement  and  snow  density  in  laboratory 
plate  sinkage  test.  (From  Shoop  and  Alger  1998.) 
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a.  At  0.097  mm/s. 


b.  At  0.98  mm/s. 


Figure  32.  Deformation  of  fresh  snow  under  a  piate.  (From  Fukue  1979.) 


Soil  model 

The  material  model  used  for  the  wheel-on- soil 
simulation  was  also  a  Drucker-Prager  cap  plasticity 
model.  Initially  a  very  soft  thawing  soil  was  to  be 
simulated,  but  the  testing  of  the  soil  parameters  be¬ 
came  quite  extensive  and  will  be  completed  in  the 
future.  A  compacted  sand  material  called  McCormick 
Ranch  sand  was  used  for  initial  modeling. 
McCormick  Ranch  sand  is  a  clayey  sand,  classified 
as  SC,  from  Albuquerque,  New  Mexico,  and  is  simi¬ 
lar  in  grain  size  distribution  to  the  Lebanon  sand  used 
in  mobility  experiments  on  frozen  and  thawing  soils 
at  CRREL  (Fig.  33).  The  McCormick  Ranch  sand 
was  tested  in  a  compacted  state,  whereas  soils  that 
have  been  subjected  to  freeze-thaw  are  usually  loose 


and  wet.  However,  the  compacted  McCormick  Ranch 
sand  model  is  useful  for  simulating  a  vehicle  on  an 
unpaved  secondary  road  under  good  (summer)  road 
conditions. 

The  mechanical  behavior  of  the  McCormick 
Ranch  sand  was  documented  by  Mazanti  and  Holland 
(1970a,  1970b),  and  material  parameters  for  a  modi¬ 
fied  Drucker-Prager  cap  model  were  developed  in 
DiMaggio  and  Sandler  (1971)  and  implemented  in 
ABAQUS  (HKS  1996).  The  material  parameters  are 
given  in  Table  4,  and  the  cap  hardening-softening 
behavior  is  defined  using  the  piecewise  linear  ap¬ 
proach  given  in  Table  5.  The  material  model  yields 
excellent  comparison  to  uniaxial  strain  test  results 
(Fig.  34). 


Grain  Size  (mm) 

Figure  33.  Grain  size  curves  for  McCormick  Ranch  sand  (Mazanti  and  Hol¬ 
land  1970a)  and  Lebanon  sand  (Shoop  1990, 1993a). 
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Table  5.  Hardening  model  for 
McCormick  Ranch  sand  (HKS 
1996). 

Hydrostatic 
stress,  Pb 
(MPa) 

Volumetric  plas¬ 
tic  strain,  8vol 

0.138 

0 

0.172 

0.005 

0.434 

0.01 

0.896 

0.02 

1.655 

0.03 

2.758 

0.04 

4.137 

0.05 

6.895 

0.06 

13.447 

0.10 

Table  4.  Material  constants  for  Drucker-Prager  cap  model  of 

McCormick  Ranch  sand  (HKS  1996). 

Material  parameter 

Value 

Young’s  modulus,  E 

689  Mpa  (100000  psi) 

Poisson’s  ratio,  v 

0.25 

Drucker-Prager  material  angle  of  friction,  p 

14.56° 

Drucker-Prager  material  cohesion,  d 

1.194  MPa  (173.2  psi) 

Cap  eccentricity  parameter,  R 

0.1 

Initial  value  of  volumetric  plastic  strain,  eJoJo 

0.001 

Flow  stress  ratio,  K 

1.0 

Axial  Strain 


Figure  34.  Validation  of  McCormick  Ranch  sand 
Drucker-Prager  cap  model  with  uniaxial  strain  test  data 
from  DiMaggio  and  Sandler  (1976).  (After  HKS  1996.) 
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TIRE  MODELS  FOR  A  DEFORMABLE 
SUBSTRATE 

Tire  model  construction 

Background 

To  apply  a  tire  model  to  deformable  terrain,  a 
model  is  needed  that  is  efficient  yet  accurately  por¬ 
trays  the  tire  structural  behavior.  Specifically,  an  ac¬ 
curate  model  of  the  contact  patch  is  critical  for  simu¬ 
lating  the  impact  of  deformable  terrain  on  tire 
performance.  Models  commonly  used  for  tire  design 
purposes  very  accurately  predict  deformation  of  the 
complete  tire,  including  the  interaction  of  the  internal 
components.  However,  since  our  concern  is  only  the 
deformation  as  it  relates  to  the  contact  region  and  its 
ability  to  roll  across  a  surface,  simpler  models  can  be 
employed  for  better  computational  efficiency.  To  this 
purpose,  several  tire  models  were  evaluated  for  effi¬ 
ciency  and  comparison  to  measured  tire  mechanical 
response,  particularly  deformation  and  stress  related 
to  contact. 

Modern  tires  are  structurally  complex,  consisting 
of  layers  of  belts,  plies,  and  bead  steel  imbedded  in 
rubber  (Fig.  7).  Materials  are  often  anisotropic,  and 
rubber  compounds  vary  through  the  tire  structure. 
Models  developed  for  tire  design  are  extremely  de¬ 
tailed,  modeling  each  material  within  the  tire  (belts, 
plies,  rubber)  for  internal  tire  stress  analysis,  wear, 
and  vibration  modal  analysis.  For  these  types  of 
models,  the  tire  may  consist  of  shell,  continuum,  and 
rebar  elements.  The  simulation  consists  of  mounting 
the  tire  on  the  rim,  inflating  it,  lowering  the  tire  onto 
a  rigid  surface,  and  applying  various  loading  condi¬ 
tions.  Such  models  are  extremely  large  and  take 
many  hours  to  run.  The  tire  properties  and  associated 
modeling  details  for  tire  models  used  by  industry  are 
extremely  difficult  to  obtain.  Among  the  more  nota¬ 
ble  published  contributions  to  tire  modeling  are  nu¬ 
merous  publications  by  Joseph  Padovan  at  the  Uni¬ 
versity  of  Akron  (Padovan  1977,  Padovan  et.  al  1992, 
Padovan  and  Padovan  1993,  1994a,  1994b)  and  mod¬ 
els  of  agricultural  tires  by  Hu  and  Abeels  (1994).  An 
example  of  such  a  model  is  the  National  Tire  Model¬ 
ing  Program  (NTMP)  radial  automobile  tire  shown  in 
Figure  35.  The  NTMP  is  a  generic  model  used  for 
calibrating  modeling  programs  across  the  tire  indus¬ 
try.  Prior  to  1997  this  type  of  model  required  over  15 
CPU  hours  on  a  SUN  SPARCstation  10  for  inflation 
and  no- slip  contact  on  a  rigid  surface  (i.e.  no  rolling). 

In  this  research  project,  four  types  of  tire  models 
were  evaluated  for  suitability  to  rolling  on  deform¬ 
able  terrain:  1)  a  rigid  tire,  2)  a  simplified  tire  model 


Figure  35.  National  Tire  Modeling 
Program  (NTMP)  model  of  a  ra¬ 
dial  tire.  (After  Danielson  and 
Noor  1997.) 

using  methodology  developed  by  Darnell  at  the  Uni¬ 
versity  of  Michigan  for  use  in  vehicle  dynamics,  3)  a 
tire  model  of  the  type  used  for  harmonic  vibration 
modal  analysis,  with  a  smooth  tread,  and  4)  a  tire 
model  similar  to  model  3  except  with  a  straight  (lon¬ 
gitudinal)  ribbed  tread.  The  modal  analysis  models 
were  provided  by  Goodyear  through  a  Cooperative 
Research  and  Development  Agreement  and  required 
only  minor  modifications  to  accommodate  surface 
contact  and  rolling.  All  of  the  models  were  built  to 
represent  the  tires  used  in  the  experimental  test  pro¬ 
gram  for  comparison  to  measured  tire  behavior  in 
terms  of  deflection,  contact  area,  deflected  sidewall 
profile,  contact  stress  distribution,  and  performance 
on  deformable  terrain  (snow  and  soil). 

Tires  simulated 

Three  tires  were  involved  in  the  data  collection 
phase  of  the  project.  The  first  tire  was  a  Michelin 
XCH4  235/75  R15.  This  tire  was  used  in  several  ex¬ 
periments  on  snow  and  soil  under  various  moisture, 
compaction,  and  temperature  conditions.  The  second 
tire,  a  Goodyear  Wrangler  AT  235/75  R15,  was  cho¬ 
sen  to  be  similar  to  the  first  in  construction  and  be¬ 
havior.  The  third  was  a  highway  tire  with  a  ribbed 
tread  pattern,  the  Goodyear  Wrangler  HT  235/75 
R15.  The  three  tires  were  similar  in  size,  construc¬ 
tion,  deflection,  and  contact  area  characteristics,  and 
data  from  all  three  were  used  in  the  following  com¬ 
parisons.  Although  the  tread  patterns  and  compounds 
varied  (Fig.  36),  which  would  strongly  influence  trac¬ 
tion,  this  study  focused  on  the  tire  deformation  and 
rolling  resistance,  so  the  tires  can  be  considered 
equivalent.  The  dimensions  of  the  test  tires  are  given 
in  Table  6. 
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Figure  36.  Tires  used  in  the  experiments:  Goodyear  Wrangier  HT, 
Micheiin  XCH4,  and  Goodyear  Wrangier  AT. 


Tabie  6.  Tire  characteristics  for  the  test  tires  at  an  average 
verticai  ioad  of  6227  N  (1400  ib).  The  XCH4  tire  data  were 
used  to  represent  the  average  vaiues. 

Inflation  pressure  contact  Deflection  Max.  deflected 
(kPa)  (psi)  area  (crrf)  (cm)  (%)  section  width  (cm) 


103 

8 

509.7 

6.3 

38 

27.9 

138 

15 

419.4 

5.4 

33 

25.9 

179 

26 

364.5 

4.6 

28 

26.0 

241 

35 

281.7 

3.1 

19 

25.1 

Unloaded  tire  diameter:  73.7  cm 
Unloaded  section  width:  22.9  cm 
Unloaded  section  height:  16.4  cm 
Tread  width:  15.7  cm 
Ratio  of  voids  to  tread  contact  area:  0.44 
Inflated  overall  section  width:  24.9  cm 
Inflated  diameter:  74.7  to  75.3  cm 


All  three  tires  were  evaluated  for  their  deforma¬ 
tion  characteristics  on  a  rigid  surface.  A  comparison 
of  the  deflection  and  contact  area  measurements  col¬ 
lected  from  various  sources  is  given  in  Figures  37 
and  38.  All  of  the  data  were  collected  after  a  tire 
break-in  period  of  operation  of  at  least  160  km.  Con¬ 
tact  area  is  calculated  from  a  static  contact  “perime¬ 
ter”  area  without  accounting  for  voids  due  to  tread 
pattern  so  that  comparisons  can  easily  be  made  be¬ 


tween  the  different  tires  and  models.  Based  on  the 
good  match  in  the  measured  tire  deformations  on  a 
rigid  surface  (deflection  and  contact  area),  the  tires 
were  assumed  to  be  sufficiently  similar  for  inter¬ 
changeable  use  in  comparing  to  model  results.  The 
normal  operating  pressure  for  the  tires  is  241  kPa  (35 
psi),  and  that  is  where  the  similarity  between  models 
and  data  is  best. 
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Figure  37.  Deflection  measurements  for  the  tires  used  in  the  experi¬ 
ments  at  inflation  pressures  of  241  kPa  (35  psi)  and  179  kPa  (26  psi). 
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Figure  38.  Comparison  of  measured  contact  area  for  the  three  test 
tires  at  241  kPa  (35  psi). 


Rigid  tire  modei 

If  the  deformable  material  is  very  soft,  as  in  fresh 
snow,  the  tire  can  be  approximated  as  a  rigid  wheel. 
The  rigid  wheel  is  modeled  as  a  rigid  analytical  sur¬ 
face  0.74  m  in  diameter  and  0.272  m  wide  (for  a  full 
wheel),  with  a  0.051 -m  radius  of  curvature  on  the  tire 
shoulder.  An  unsprung  mass  of  636  kg  is  placed  at 
the  axle  to  simulate  the  weight  of  the  tire  and  the  load 
of  the  vehicle  on  the  tire  (6240  N);  the  rotational 
inertia  of  the  wheel  was  set  to  2.15  kg  m^. 

Darneii  tire  modei 

To  adequately  model  the  dynamics  of  a  rolling 
tire,  a  very  efficient  tire  model  is  needed.  Thus,  the 
tire  model  developed  by  Darnell  and  Hulbert  at  the 
University  of  Michigan  is  a  promising  alternative. 
The  Darnell  tire  model  (Darnell  et  al.  1997,  Darnell, 
in  progress)  is  an  extension  of  a  two-dimensional  tire 
modeling  concept  developed  by  Mousseau  (Mous- 
seau  and  Hulbert  1996)  for  interfacing  with  vehicle 
dynamics  simulations.  The  tire  is  modeled  using  ma¬ 
terial  properties  for  the  composite  material  rather 
than  the  individual  components.  The  composite  prop¬ 
erties  were  obtained  by  sectioning  the  tire  and  per¬ 
forming  laboratory  tests  on  each  of  the  major  section 
components  (tread,  sidewall,  and  shoulder).  A  similar 
approach  was  also  developed  at  the  University  of 
Birmingham  in  the  UK  for  tire  modal  analysis  (Burke 
and  Olatunbosun  1997a,  1997b,  1997c).  Additional 
work  on  the  Darnell  model  has  been  toward  paramet¬ 


ric  analysis  and  refinement  of  the  code  (Alverez  Sanz 
1999). 

The  objectives  of  the  Darnell  tire  model  are  simi¬ 
lar  to  those  needed  for  this  study:  an  efficient,  accu¬ 
rate,  three-dimensional  model  to  predict  spindle 
(axle)  forces  but  not  the  internal  stresses  in  the  tire. 
The  contact  forces  are  calculated  in  the  process,  so 
the  model  is  also  suitable  for  efficient  modeling  of 
the  tire-terrain  interaction,  where  accurate  modeling 
of  the  contact  stress  distribution  is  imperative. 

The  tire  is  simulated  with  three  modeling  ele¬ 
ments:  the  tread  elements  are  six-degree-of-freedom 
shell  elements  (S4R);  the  sidewall  elements  are  three- 
node,  user-defined  elements  based  on  pre-calculated 
results  from  a  sidewall  model;  and  the  contact  is 
modeled  using  a  “hard”  contact,  penalty  method 
(HKS  1998).  The  roadway  is  modeled  as  a  rigid  sur¬ 
face. 

The  sidewall  model  is  a  special-purpose  finite 
element  code  that  generates  a  lookup  table  of  side- 
wall  forces  (and  geometry)  using  ten  equal-length, 
non-extensible  beam  elements  with  inflation  pressure 
as  a  follower  load.  The  sidewall  model  is  generated 
from  the  geometry  of  the  uninflated  sidewall,  the 
bending  and  shear  stiffness  for  each  element,  and  the 
tire  inflation  pressure.  The  results  of  the  sidewall 
model  are  implemented  in  the  tire  model  through  the 
user-defined  element,  which  attaches  to  the  tread 
section  with  rotational  springs  (Fig.  39). 
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d2 


Figure  39.  Construction  eiements  of  the  Darneii  tire  modei.  (After  Darneii  et  ai.  1997.) 


The  Darnell  tire  model  was  validated  for  a  Ford 
Taurus  tire  for  lateral  forces  and  deflection  on  a 
smooth  surface  and  while  rolling  over  an  angled  ver¬ 
tical  step  (Darnell  et  al.  1997)  and  for  slip  angle  and 
aligning  moment.*  The  model  results  sufficiently 
mimic  tire  behavior  while  being  computationally 
efficient,  indicating  its  potential  for  simulations  of  a 
tire  on  deformable  terrain.  Some  alteration  in  the 
sidewall  model  may  be  necessary  to  account  for  the 
interaction  of  the  sidewall  with  the  deformable  terrain 
if  sinkage  is  substantial. 

The  same  tire  used  to  gather  motion  resistance 
data  on  various  terrain  (snow,  freezing  and  thawing 
ground)  was  cut  into  four  sections:  sidewall,  lateral 
tread,  longitudinal  tread,  and  shoulder  (Fig.  40).  The 
sidewall  was  further  dissected  into  smaller  sections 
for  future  work  expanding  the  sidewall  model.  The 
sections  were  tested  in  the  laboratory  to  determine 
the  beam  stiffness  {ET)  by  clamping  an  end  of  the  test 
section  onto  a  workbench  to  form  a  cantilever  beam. 
Deflections  were  measured  for  a  range  of  loads,  and 
stiffness  was  calculated  using  the  equation  for  a  can¬ 
tilever  beam  with  a  point  load  on  the  free  end: 

-  1 

ymsLX  ^ 

3  El 

where  ymax  =  maximum  deflection  at  the  free  end 
of  the  beam 
W  =  applied  load 
E  =  Young’s  modulus  of  the  beam 


*  Personal  communication  with  I.  Darnell,  University 
of  Michigan,  1997. 


Figure  40.  Cut  tire  sections  used  to  determine 
composite  tire  properties. 


1  =  moment  of  inertia 
I  =  beam  length. 

Young’s  modulus  {E)  is  calculated  from  El  using 
1  =  ’^beam^Lm  ^12,  whcrc  Wbeam  and  ffeeam  arc  the  width 
and  thickness  of  the  tire  samples,  respectively.  The 
test  results  are  summarized  in  Table  7. 

The  sidewall  was  divided  into  ten  equal-length 
beam  elements  in  a  special-purpose  finite  element 
code.  This  code  generates  the  sidewall  force  and  ge¬ 
ometry  table  for  the  user-defined  element  represent¬ 
ing  the  sidewall  in  the  tire  model.  The  stiffness  and 
curvature  of  the  sidewall  model  are  given  in  Table  8 
and  illustrated  in  Figure  41. 

Additional  measurements  needed  for  the  model 
are  sidewall  curvature,  tread  hoop  stiffness,  rotational 
spring  constants,  reference  angles,  tire  and  rim  geo- 
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Table  7.  Bending  stiffness  of  tire  sections. 


Width 

(cm) 

Length 

(cm) 

Thickness 

(cm) 

El 

(Nnf) 

Sidewall  section  A1 

5.2 

3.3 

1.3 

0.1564 

Sidewall  section  A2 

5.2 

2.0 

1.3 

0.0287 

Sidewall  section  A1+A2 

5.2 

5.3 

1.3 

0.1521 

Sidewall  section  B 

5.8 

6.1 

0.6 

0.0459 

Sidewall  section  C 

7.0 

2.1 

1.1 

0.0172 

Lateral  tread  section 

7.0 

16.5* 

1.7 

(0.9  w/o  tread) 

0.1550 

Longitudinal  tread  section 

Shoulder 

12.5 

27.2 

33.5  cm 
along  tread 

1.7 

(0.9  w/o  tread) 

0.5137 

*  23.4  cm  including  shoulders 


Table  8.  Properties  of  sidewall  elements  and 
tire  model  components. 

Curvature  Stiffness  per  unit 

Eiement  ( 1/m)  width  (N  m) 


1  (Bead  end) 

0.01 

2  (rigid) 

2 

6.76 

1.5 

3 

7.63 

0.15 

4 

9.95 

0.04 

5 

8.62 

0.04 

6 

9.71 

0.04 

7 

9.01 

0.04 

8 

9.09 

0.15 

9 

14.49 

0.5 

10  (Tread  end) 

17.54 

1.0 

Total  length  of  sidewall  =  16.1  cm,  each  element  L=  1.6 
cm 

Tread  attachment  angle  =  133° 

Rim  attachment  angle  =  23° 

Tread  diameter  to  end  of  last  sidewall  element  =  72.5  cm 
Tread  diameter  to  center  of  belts  =  69.5  cm 
Rim  diameter  to  first  sidewall  element  =  41 .6  cm 
Tread  longitudinal  shell  flexural  rigidity  for  unit  section 
width,  £i/  =  4.1261  Nm^ 

For  1 .1 1-mm  shell  thickness,  E,  =  3.62  x  10'°  N/m^ 
Tread  lateral  shell  flexural  rigidity  for  a  unit  section  width, 
£2/  =  2.259  Nm^ 

For  1.11  -mm  shell  thickness,  E2  =  1 .98  x  1 N/m^ 


metry,  and  sidewall  and  shoulder  stiffness.  These  are 
listed  at  the  bottom  of  Table  8. 

For  the  elements  representing  the  tread,  the  shell 
stiffness  was  calculated  with  Ei  being  along  the  tread 
in  the  longitudinal  direction  and  E2  being  across  the 
tread  in  the  lateral  direction.  An  equivalent  shell 
thickness  of  1.11  mm  was  used  in  the  model  since 
this  was  found  to  yield  a  reasonable  hoop  stiffness 
based  on  measured  changes  in  circumference  with 
inflation  pressure.* 


*  Personal  communication  with  I.  Darnell,  University 
of  Michigan,  2000. 


Figure  41.  Tire  geometric  parameters  for 
the  sidewaii  modei. 
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Viewport:  1  Model:  Model-1  Module:  Visualization 


Viewport:!  Model: Model-1  Module: Visualization 


a.  Coarse  (3  x  36)  mesh. 

Figure  42.  Shoop-Darnell  tire  modei  (user  eiements  not  shown). 


The  Shoop-Darnell  model  (the  model  of  the  test 
tire  created  with  the  Darnell  methodology)  was  con¬ 
structed  using  ABAQUS  for  a  quasi-static  simulation. 
This  adequately  simulates  low- speed  vehicle  travel 
and  is  appropriate  for  5  mph,  which  is  the  standard 
test  speed  for  off-road  traction  and  motion  resistance 
testing.  The  tire  model  is  first  subjected  to  internal 
pressure  (inflated),  then  lowered  onto  the  surface  and 
allowed  to  equilibrate  with  the  applied  vertical  load 
before  rolling.  The  Darnell  methodology  for  the  tire 
model  was  validated  for  the  test  tire  on  a  hard  surface 
and  then  used  to  simulate  the  three-dimensional  con¬ 
tact  between  two  deformable  bodies  on  a  homogene¬ 
ous  terrain,  first  using  an  elastic  terrain  model  and 


then  a  plasticity  model  for  the  terrain  material.  Two 
Shoop-Darnell  mesh  configurations  were  used  in  the 
analysis:  a  coarse  and  a  fine  mesh  (Fig.  42).  Since  the 
sidewalls  are  modeled  with  user-defined  elements, 
they  do  not  appear  on  the  standard  finite  element 
visualization. 

Modal  analysis  tire  model:  Smooth  tread 

The  smooth  tread  tire  model  is  typical  of  the  type 
of  models  used  for  tire  vibration  modal  analysis.  The 
carcass  is  composed  of  a  single  layer  of  four-node, 
reduced-integration  shell  elements  (S4R)  with  mate¬ 
rial  properties  representing  the  composite  behavior 
through  the  carcass  thickness.  The  treadcap  is  con- 
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Figure  42  (cont.).  Shoop-Darnell  tire  modei. 


stmcted  of  linear,  hybrid  continuum  elements 
(C3D8H),  with  constant  pressure  (simulating  the 
nearly  incompressible  nature  of  rubber).  The  result  is 
a  reasonable  approximation  of  both  the  structural 
behavior  and  the  contact  patch.  The  model  is  con¬ 
structed  using  the  cross  section  shown  in  Figure  43, 
repeated  every  three  degrees  around  the  axis  of  the 
hub  (creating  120  cross  sections).  For  most  of  the 
simulations  the  tire  model  was  cut  in  half  using  a 
vertical  symmetry  plane  along  the  longitudinal  axis. 
The  wheel  mass  is  represented  by  one  mass  element 
and  one  rotary  inertia  element  located  at  the  wheel’s 
centerpoint  on  the  symmetry  plane. 


The  analysis  is  done  in  three  steps:  inflating  the 
tire  to  the  desired  pressure,  vertically  loading  it  based 
on  the  vehicle’s  weight,  and  rolling  it  along  the  road 
by  translating  its  centerpoint  while  allowing  free  rota¬ 
tion  about  its  axis.  The  creation  and  seating  of  the 
bead  are  accomplished  by  multi-point  constraint  and 
fixed  displacement  of  the  first  four  rings  of  nodes 
closest  to  the  wheel  rim.  During  the  inflation  and 
loading  steps,  no  friction  is  included  on  the  contact 
surface.  Variations  of  the  analysis  include  a  static 
rolling  step,  a  dynamic  rolling  step,  and  changes  in 
the  definition  of  the  shell  element  local  coordinate 
system. 
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Viewports  Model:  Model-1  Module:  Visualization 
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P235/75R1S  wrangler  KT, 

ODB;  1490_30.o<U)  ABAQUS/Standard  5.8-14  Fri  Apr  14  16:26:25  EDT  2000 
Step:  Step  3  Zncrement  20:  Step  Time  =  1.000 


Figure  43.  Modal  analysis  tire  model  with  smooth  tread. 
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Modal  analysis  tire  model:  Straight  ribbed  tread 

The  ribbed  tread  tire  model  is  also  a  modal  analy¬ 
sis  tire  model.  Construction  is  similar  to  the  smooth 
tread  model,  with  the  carcass  composed  of  four-node 
reduced-integration  shell  elements  and  the  treadcap 
constructed  of  linear,  hybrid  continuum  elements 
with  constant  pressure.  The  major  difference  is  the 
four  grooves  across  the  width  of  the  tire  created  by 
omitting  sections  of  the  tread  continuum  elements 
(Fig.  44).  The  wheel  is  represented  by  one  mass  ele¬ 
ment  and  one  rotary  inertia  element  located  at  the 
centerpoint.  Details  for  the  smooth  and  ribbed  tread 
models  are  given  in  Table  9. 

The  analysis  is  completed  in  eleven  steps,  com¬ 
prising  inflation,  bead  seating,  lowering  to  the  road 
surface,  and  vertical  loading.  The  loading  occurs  in 
eight  steps  of  alternating  loose  and  tight  convergence 
control.  Inflation  is  handled  similarly.  The  creation 
and  seating  of  the  bead  are  accomplished  using  multi¬ 
point  constraints  on  the  first  five  rings  of  nodes  clos¬ 
est  to  the  wheel  rim.  An  additional  step  translates  the 


wheel’s  centerpoint  while  allowing  the  tire  to  freely 
rotate  about  its  axis  due  to  friction,  simulating  a 
towed  wheel.  The  contact  is  frictionless  during  the 
loading  process,  and  friction  is  added  during  the  roll¬ 
ing  step.  Variations  of  the  analysis  include  a  static  or 
dynamic  analysis  procedure  during  rolling  and 
changes  in  the  convergence  controls. 

The  towing  phase  of  the  dynamic  analysis  is  bro¬ 
ken  into  two  steps.  The  first  is  to  accelerate  the  tire 
from  rest  to  1  m/s  at  a  constant  acceleration  of  1  m/s^ 
for  one  second  (traversing  500  mm).  The  second  step 
continues  to  roll  the  tire  at  a  constant  velocity  of  1 
m/s  for  two  more  seconds  (2000  mm).  For  all  cases 
the  analyses  converge  for  a  while  but  eventually 
reach  a  point  where  the  tire  stops  moving  by  taking 
increasingly  smaller  time  steps  due  to  convergence 
difficulties.  The  problem  occurs  when  the  shell  is 
subjected  to  simultaneous  distortion  and  rotation  as 
the  rotation  angle  approaches  360°.  This  is  a  bug  in 
the  ABAQUS  shell  element  that  will  be  corrected  in 
future  versions  (Bug  Report  v58_2449). 


Figure  44.  Modal  analysis  tire  model  with  ribbed  tread. 
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Figure  44  (cont.)-  Modal  analysis  tire  model  with  ribbed  tread. 


Table  9.  Features  of  the  modal  analysis  models  with  smooth  and  straight 
ribbed  tread. 


Smooth  tread 

Half  carcass  model 

Straight  ribbed  tread 
Half  carcass  model 

Carcass  elements  (number  and  type) 

3360  S4R 

6240  S4R 

Treadcap  elements  (number  and  type) 

1440  C3D8H 

2400  C3D8H 

120  C3D6H 

120  C3D6H 

Tread  rubber:  neo-Hookean,  Cio 

1.109596 

2.408340 

Tread  density  (Ns^/mm"^) 

0.814575E-9 

6.6255E-10 

Tread  gaps  (width  x  depth,  mm) 

None 

5.1  X  16.6 

5.3  X  16.2 

Total  no.  elements 

6602 

13562 

Total  model  variables 

32172 

63396 
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Evaluation  of  tire  models  on  a  rigid  surface 

The  tire  models  were  applied  to  a  rigid  surface 
through  a  range  of  pressures  and  loads  for  compari¬ 
son  to  measured  deflections,  contact  patch  areas,  di¬ 
mensions,  and  stresses. 

Deflection 

A  complete  data  set  for  deflection  was  gathered 
for  the  Wrangler  AT  (Fig.  45)  for  use  in  comparisons 
with  model  results  from  the  Shoop-Darnell  model, 
the  smooth  tread  model,  and  the  ribbed  tread  model. 
Each  model  was  loaded  at  a  range  of  vertical  loads 
(from  0  to  8000  N)  at  three  tire  inflation  pressures: 


241  kPa  (35  psi),  179  kPa  (26  psi),  and  103  kPa  (15 
psi).  The  suggested  inflation  pressure  for  these  tires  is 
241  kPa.  Lower  inflation  pressures  are  sometimes 
used  when  driving  in  off-road  and  marginal  road 
conditions  (snow  and  soft  soil)  and  for  minimizing 
damage  to  unpaved  travel  surfaces.  Therefore,  the 
two  lower  inflation  pressures  were  also  evaluated. 
Performance  at  a  range  of  inflation  pressures  is  also 
of  interest  to  industries  using  vehicles  with  Central 
Tire  Inflations  Systems  (CTIS)  (i.e.  military,  forestry, 
and  agriculture). 

A  compilation  of  the  data  and  model  results  for 
each  inflation  pressure  is  shown  in  Figure  46.  The 


Figure  46.  Comparison  of  measured  and  modeied  defiection  for  the 
Shoop-Darneii,  smooth  tread,  and  ribbed  tread  modeis  for  three  infia- 
tion  pressures. 
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Vertical  Load  (N) 

Figure  47.  Comparison  of  aii  measured  and  modeied  defiections  at 
241  kPa  (35  psi)  infiation  pressure. 


Figure  48.  Measured  contact  areas  at  three  infiation  pressures. 


model  predictions  follow  the  same  general  trend  as 
the  data  for  all  tire  pressures,  with  the  largest  differ¬ 
ence  occurring  at  the  lower  inflation  pressure.  This  is 
understandable  since  undesirable  tread  and  sidewall 
behavior  can  occur  when  tires  are  underinflated,  and 
these  pressures  and  behaviors  are  neither  within  the 
design  range  nor  accurately  accounted  for  in  the 
models.  On  the  other  hand,  the  agreement  between 
data  and  models  at  the  standard  inflation  pressure  is 
exceptionally  good,  as  illustrated  in  more  detail  in 
Figure  47,  showing  all  measured  and  modeled  deflec¬ 
tions  at  241  kPa. 

Contact  area 

Comparisons  of  contact  area  were  done  much  the 
same  as  the  deflection  comparisons.  The  measured 


contact  areas  for  the  three  inflation  pressures  are 
given  in  Figure  48.  All  contact  areas  are  based  on  the 
perimeter  of  the  contact,  without  accounting  for  voids 
within  the  area  due  to  tread  design.  Comparisons  of 
the  model  results  (the  Shoop-Darnell  model  with 
both  coarse  and  fine  mesh,  the  smooth  tread  model, 
and  the  ribbed  tread  model)  with  measured  data  are 
shown  in  Figure  49. 

In  general,  the  agreement  between  the  models  and 
data  is  reasonable.  Surprisingly  the  Darnell  gives  the 
best  fit,  but  the  model  mesh  is  coarse,  which  makes 
the  results  erratic  since  the  contact  area  changes 
abruptly  as  each  shell  element  makes  contact.  The 
refined  mesh  on  the  Darnell  model  is  only  a  slight 
improvement,  and  the  results  are  worse  at  the  higher 
loads  where  the  contact  nodes  pop  off  the  surface  just 
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Figure  49.  Comparison  of  measured  and  modeied  contact  areas  at  three 
infiation  pressures. 
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Figure  50.  Comparison  of  Shoop-Darneii  modei  resuits  with  measured 
contact  areas  (for  coarse  and  fine  meshes). 
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Figure  51.  Comparison  of  the  smooth  and  ribbed  modai  anaiysis  tire  modeis 
with  measured  contact  areas. 

inside  of  the  tire  shoulder  due  to  buckling  of  the  car¬ 
cass  (as  discussed  later).  This  behavior  strongly  af¬ 
fects  the  contact  stress  and  area  since  adjacent  nodes 
may  alternate  between  high  loads  and  no  contact. 

Figure  50  compares  the  Shoop-Darnell  models  to  the 
measured  data. 

Comparisons  of  the  smooth  and  ribbed  tread  mo¬ 
dal  analysis  models  to  measured  data  are  shown  in 
Figure  51.  Both  models  provide  contact  areas  less 
than  the  measured  data  at  the  high  loads,  indicating 
that  the  models  are  slightly  stiff.  The  ribbed  model 


contact  area  is  slightly  less  than  the  smooth  model, 
even  after  being  corrected  for  the  tread  voids. 

Contact  stress  distribution 

Contact  stress  measurements  for  the  Wrangler  HT 
were  provided  by  Goodyear.  For  comparison,  the 
modeled  contact  stresses  were  displayed  at  the  same 
contour  levels  and  color  scheme,  as  shown  in  Figure 
52.  Both  models  and  measurements  show  high  stress 
values  at  the  tire  shoulder  and,  to  a  lesser  extent, 
along  the  tire  centerline.  Agreement  with  the  meas¬ 
ured  contact  stress  distribution  is  excellent. 
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Figure  52.  Measured  and  modeled  contact  stress  distribution  for  half  carcass  on  a  hard  surface 
(207  kPa  inflation  and  6627  N  load). 


Figure  53.  Irregular  stress  contours  generated  by 
the  Shoop-Darnell  model  on  a  rigid  surface  (full 
carcass  width). 


Figure  54.  Close-up  view  of  the  contact  for  the 
ribbed  tire  model,  showing  buckling  just  inside 
the  shoulder  area  (half  tire  carcass). 


During  contact  on  a  rigid  surface,  the  Shoop- 
Darnell  model  suffered  from  buckling  just  inside  the 
tire  shoulder  area,  which  strongly  affected  the  contact 
stress.  Although  a  finer  mesh  would  normally  smooth 
the  contact  stress  contours,  the  buckling  causes  alter¬ 
nate  nodes  to  lift  off  the  surface,  resulting  in  the 
stress  pattern  shown  in  Figure  53.  Eventually,  pro¬ 


gressively  finer  meshes  fail  to  converge  [similar  be¬ 
havior  was  noted  and  explored  by  Alvarez  Sanz 
(1999)].  The  buckling  behavior  also  occurs  to  a  lim¬ 
ited  extent  in  the  physical  tire  and  in  the  modal 
analysis  tire  models  (seen  in  the  contours  in  Figure 
52).  The  buckling  is  worse  at  low  inflation  pressures 
and  high  loads  (Fig.  54).  The  impact  is  minor  for  the 
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Figure  55.  Measured  and  modeled  sidewall  profiles  for  the  Shoop- 
Darnell  tire  model. 


actual  tires  and  the  models  where  the  rubber  tread 
elements  conform  to  the  surface,  softening  the  con¬ 
tact.  For  our  purposes,  this  behavior  disappears  when 
the  tire  model  is  applied  to  a  terrain  surface,  which 
softens  the  contact  through  terrain  compliance. 

Sidewall  profile 

Measurements  of  the  deflected  shape  of  the  side- 
wall  under  various  loads  and  inflation  pressures  were 
collected  to  check  the  calibration  of  the  Shoop- 
Darnell  tire  sidewall  model.  Selected  measurements 
of  the  profiles,  along  with  comparisons  to  the  mod¬ 
eled  geometry,  are  given  in  Figure  55.  The  modeled 
sidewall  geometry  closely  duplicates  the  measured 
geometry.  The  modeled  sidewall  is  slightly  longer 
than  the  measured  profiles  since  the  measurements 
did  not  include  the  sidewall  portion  covered  by  the 
wheel  rim. 

Hard-surface  rolling  resistance 

Ideally  forces  generated  from  rolling  the  model 
should  be  comparable  to  the  measured  hard-surface 
rolling  resistance.  However,  rolling  resistance  is 
largely  due  to  viscoelastic  behavior  in  the  rubber 
compounds,  and  these  models  omit  the  viscoelastic 
material  properties  for  increased  efficiency.  There¬ 
fore,  the  modeled  hard-surface  rolling  resistance 
forces  reflect  only  resistance  due  to  interface  shear. 
For  an  interface  friction  coefficient  of  0.825,  similar 
to  asphalt  pavement,  the  rolling  forces  modeled  are 
less  than  10  N  at  an  inflation  pressure  of  241  kPa, 
while  measured  values  are  on  the  order  of  25  N  [at 
241  kPa  and  low  speed  (8  kph)]. 


Run  times 

Model  efficiency  is  of  primary  importance  for  use 
in  rolling  a  tire  on  a  deformable  material,  so  the  mod¬ 
els  were  run  with  various  options  and  using  different 
computers  in  order  to  assess  the  effects  on  run  time. 
The  primary  variables  considered  were  the  effect  of 
the  tread  (smooth  vs.  ribbed),  static  or  dynamic 
analysis,  convergence  controls,  optimum  use  of  sym¬ 
metry,  and  computing  system.  The  results  of  the 
many  runs  are  included  in  the  following  discussions 
and  qualitatively  summarized  in  Table  10.  Examples 
of  some  of  the  run  times  are  given  at  the  bottom  of 
the  table.  Consideration  was  also  given  to  how  far  the 
tire  would  roll  in  the  model  run.  A  roll  distance  of 
between  one  third  and  two  revolutions,  depending  on 
the  tire,  is  needed  to  remove  abnormal  stresses  built 
up  in  the  tire  structure  due  to  lowering  it  onto  the 
surface.  One  revolution  was  considered  the  minimum 
requirement.  All  comparison  simulations  were  done 
using  the  implicit  solver. 

Static  dynamic.  The  static  models  were  run  in  a 
steady-state  type  of  analysis.  The  same  models  were 
modified  for  a  dynamic  analysis  where  the  inertia  of 
the  system  is  considered.  In  general,  the  static  models 
roll  over  twenty  times  farther  in  less  than  twice  the 
CPU  time  than  the  dynamic  models.  However,  the 
dynamic  models  used  less  memory.  Ultimately,  a 
dynamic  model  is  the  desired  choice  to  include  mass, 
velocity,  and  acceleration  effects  and  for  operating 
the  combined  model  using  an  explicit  finite  element 
code. 

Full  tire  half-tire.  The  full-tire  version  of  the 
ribbed  model  is  prohibitively  expensive  in  terms  of 
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Table  10.  Qualitative  summary  of  runtime  analysis. 


Model  parameters  compared 


Change  in  performance 


Static  vs.  dynamic  analysis 

Half  tire  carcass  vs.  full  tire  carcass 
Smooth  vs.  ribbed  tread  model 
SGI  Octane  vs.  Cray  J932 
Tight  vs.  loose  convergence  controls 


Static  models  rolled  2  m  vs.  0.08  m  before  conver¬ 
gence  problems 
Half  tire  runs  three  times  faster 
Smooth  model  runs  four  times  faster 
Octane  is  two  to  three  times  faster 
Loose  controls  run  faster  but  tight  controls  let  the  tire 
roll  further  before  convergence  problems 


Example  run  times 

Ribbed  tread,  full  tire,  tight  convergence  Rolls  25  mm  in  488  hr  on  J932 

Ribbed  tread,  half  tire,  dynamic  analysis  Rolls  46  mm  in  252  hr  on  J932 

Smooth  tread,  half  tire,  static  analysis  Rolls  2  m  in  22  hr  on  Octane 

Shoop-Darnell  model  (full  tire) Roll  20  m  in  11  min  on  Octane 


run  time.  It  takes  almost  three  times  as  long  as  the 
half-tire  model  to  go  the  same  distance.  No  full-tire 
model  has  stopped  due  to  non-convergence,  but 
rather  due  to  system  errors  or  time  limits.  In  the  fu¬ 
ture,  if  side  forces,  cornering  forces,  or  steering  reac¬ 
tions  are  investigated,  then  a  full-tire  model  will  be 
required. 

Smooth  tread  ribbed  tread.  The  model  with  the 
smooth  tread  has  fewer  nodes  and  elements,  is  less 
geometrically  complex,  and  uses  default  convergence 
controls.  It  runs  more  than  four  times  faster  than  the 
ribbed  tread  model  during  the  static  analysis  and 
more  than  eleven  times  faster  for  the  dynamic  analy¬ 
sis.  In  most  instances  the  accuracy  of  the  smooth 
model  is  as  good  as  or  better  than  the  ribbed  model 
(i.e.  the  model  results  are  closer  to  the  experimental 
data).  In  the  dynamic  analyses  the  smooth  tread 
model  went  almost  twice  as  far  as  the  ribbed  tread 
model. 

Cray  J932  SGI  octane.  The  models  are  being 
run  on  an  SGI  Octane  (Intrepid,  120  megaflop)  and  a 
Cray  J932  (Chilkoot,  12  CPU  at  2.4  gigaflop  each, 
200  total  if  optimized  for  this  platform).  Generally 
the  same  models  took  two  to  three  times  longer  to  run 
on  the  Cray  J932  than  on  the  SGI  Octane.  This  is 
likely  because  HKS  (the  ABAQUS  software  devel¬ 
oper)  uses  the  SGI  as  a  development  platform,  which 
results  in  increased  efficiency  on  that  equipment. 

Convergence  controls.  For  the  static  analysis 
ribbed  tire  model,  the  tire  is  maneuvered  in  several 
steps,  with  the  convergence  controls  adjusted  at  each 
step  to  increase  efficiency.  The  loose  convergence 
controls  will  allow  the  model  to  run  faster,  but  con¬ 
vergence  stops  after  rolling  only  a  short  distance.  For 
dynamic  rolling  the  model  using  the  default  conver¬ 
gence  controls  runs  faster,  and  the  tire  rolls  three 
times  farther,  than  when  convergences  controls  are 
looser  or  tighter  than  default  values. 


COMBINED  TIRE-TERRAIN  MODEL 
Modeling  the  tire-terrain  contact  interface 

Contact 

To  merge  the  tire  and  terrain  models,  the  interface 
between  the  two  meshes  must  interact.  This  is  ac¬ 
complished  by  defining  where  contact  is  allowed, 
how  the  contact  occurs,  and  how  forces  are  transmit¬ 
ted.  For  two  meshes  to  come  into  contact,  a  surface 
must  be  placed  on  each  material.  For  a  deformable 
tire  on  snow  or  soil,  one  surface  covers  the  outside  of 
the  tire  and  another  surface  lies  on  the  top  of  the  ter¬ 
rain  layer.  The  contact  is  defined  by  the  proximity  of 
the  surfaces.  When  both  material  meshes  are  deform¬ 
able,  one  surface  is  chosen  as  a  “master”  surface  and 
the  other  as  a  “slave”  surface.  The  more  finely 
meshed  surface  is  usually  chosen  as  the  slave  since 
the  nodes  on  the  slave  surface  are  not  allowed  to 
penetrate  the  master  surface. 

Contact  can  be  enforced  using  either  “penalty”  or 
“kinematic”  contact  methods.  Both  methods  were 
used.  Penalty  contact  was  implemented  for  pseudo¬ 
static  analyses,  where  inertia  forces  were  not  consid¬ 
ered.  In  penalty  contact  the  penetration  of  the  sur¬ 
faces  into  each  other  is  resisted  by  linear  spring 
forces  with  values  proportional  to  the  distance  of 
penetration.  These  forces  pull  the  surfaces  into  an 
equilibrium  position  with  no  penetration.  Kinematic 
contact  was  used  in  some  of  the  dynamic  simulations 
(using  the  explicit  code).  Kinematic  contact  considers 
the  inertia  forces  of  the  material  when  calculating  the 
forces  relative  to  the  surface  positions.  Although  the 
overall  model  results  were  the  same  with  both  contact 
methods,  the  kinematic  contact  was  slightly  more 
representative  of  field  observations  of  the  snow  de¬ 
formation  beneath  the  wheel  (Haehnel  2000). 
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Friction 

Forces  parallel  to  the  interface  are  transmitted 
based  on  a  friction  law.  Friction  is  a  complex  phe¬ 
nomenon  involving  stress  and  strain  distributions; 
heat  transfer;  hydrodynamic  and  elastohydrodynamic 
fluid  flow;  material  structure;  chemical  interactions 
between  surfaces,  surface  coatings,  and  lubricants; 
and  phase  change.  Current  models  of  friction  are  em¬ 
pirical  or  semi-empirical  formulations  in  which  sev¬ 
eral  variables  and  combinations  of  variables  are 
modified  with  coefficients  and  exponents  based  on 
experimental  observations.  These  equations  apply 
only  to  the  system  for  which  they  were  developed, 
though  they  may  be  used  to  estimate  the  behavior  of 
similar  systems.  They  are  not  applicable  to  other  ex¬ 
perimental  configurations  or  test  conditions  and  can¬ 
not  be  extrapolated  outside  of  their  bounds  with  any 
degree  of  confidence.  Thus,  an  accurate  model  of 
friction  is  currently  non-existent  and  is  not  likely  to 
be  developed  in  the  immediate  future.  The  only  rea¬ 
sonable  way  to  predict  the  behavior  of  most  frictional 
systems  today  is  to  test  them  (Ludema  1996a,  1996b, 
Barber  1991). 

The  experimental  measurement  of  tire-terrain  fric¬ 
tion  (including  pavement)  is  typically  performed  as  a 
vehicle  traction  or  braking  test  [a  discussion  of  winter 
traction  test  methods  is  available  in  Shoop  et  al. 
(1994)].  Test  data  are  generally  reported  as  a  traction 
(or  friction)  coefficient  (a  peak  value  or  an  averaged 
value  over  a  specified  range)  or  as  a  traction  curve 
with  respect  to  wheel  slip.  Traction  coefficients  and 
curve  shapes  are  functions  of  both  the  tire  and  the 
surface.  Curves  of  traction  data  with  wheel  slip  can 
be  implemented  in  ABAQUS  or  other  finite  element 
analysis  codes  by  specifying  friction  as  a  function  of 
relative  velocity.  Liu  and  Wong  (1996)  took  this  ap¬ 
proach  using  a  friction  law  of  the  following  form  and 
implementing  it  in  a  tire-terrain  model  using  the  gen¬ 
eral-purpose  finite  element  program  called  MARC: 

T  =  pa(l-e“-’/^)  (26) 

where  jl  =  friction  coefficient 
X  =  shear  stress 
a  =  normal  stress 

j  =  relative  slip  distance  between  the 
wheel  and  the  terrain 

k  =  constant  having  the  same  units  as  j. 

This  equation  is  of  the  type  commonly  used  in 
conjunction  with  soil  shear  data.  However,  Liu  and 
Wong  were  not  convinced  that  this  equation  ade¬ 
quately  described  traction. 


The  simulation  of  driving  and  braking  traction 
would  be  an  important  extension  of  this  research. 
More-elaborate  friction  models  could  be  generated 
based  on  the  multitude  of  experimental  traction  data 
obtained  by  CRREL.  However,  this  need  for  a  trac¬ 
tion  curve  a  priori  seems  to  negate  the  usefulness  of 
the  model  for  predicting  traction  (i.e.  the  traction 
input  is  needed  to  get  a  valid  traction  output).  Thus, 
the  current  modeling  effort  focused  on  the  tire  forces 
developed  due  to  deformation  of  the  snow  rather  than 
due  to  the  interaction  at  the  interface  (i.e.  a  model  of 
rolling  resistance  rather  than  traction).  Therefore,  for 
this  project  the  simple  Coulomb  friction  model  was 
used: 

l^  =  %-  (27) 

Rigid  wheei  in  snow 

Model  construction 

The  snow-wheel  model  simulates  one  wheel  of 
the  CRREL  Instrumented  Vehicle  moving  through  a 
range  of  shallow  (20  cm)  to  deep  (50  cm)  fresh  snow 
with  a  density  of  200  kg/m  .  The  wheel  is  modeled  as 
a  rigid  surface  0.74  m  in  diameter  and  0.272  m  wide 
(for  a  full  wheel)  with  a  0.05 1-m  radius  of  curvature 
on  the  tire  shoulder.  An  unsprung  mass  of  636  kg  is 
placed  at  the  hub  to  simulate  the  weight  of  the  vehi¬ 
cle  and  tire;  the  rotational  inertia  of  the  wheel  was  set 
to  2.15  kg  m^.  The  snow  was  modeled  with  eight- 
node  continuum  elements  near  the  wheel  and  four- 
node  infinite  elements  for  the  far-field  snow  cover. 
The  Drucker-Prager  cap  material  model  with  linear 
elasticity  was  used  for  the  snow,  and  a  linear  elastic 
model  was  used  for  the  infinite  elements.  Coulomb 
friction  was  applied  at  the  tire-snow  interface  (jl  = 
0.3).  The  deep  snow  was  modeled  using  Adaptive 
Lagrangian-Eularian  (ALE)  meshing  to  accommo¬ 
date  the  large  displacements,  so  slide  planes  were  not 
needed  along  the  sides  of  the  tire.  The  shallow  snow 
model  had  11,520  continuum  elements  and  1,418 
infinite  elements;  the  deep  snow  model  had  14,400 
continuum  and  1,962  infinite  elements.  In  all  models 
the  snow  rested  on  a  rigid  surface.  The  shallow  snow 
model  is  shown  in  Figure  56. 

The  rolling  resistance  test  is  simulated  by  first 
lowering  the  wheel  into  the  snow  by  gravity,  accelerat¬ 
ing  it  to  the  desired  speed,  and  then  translating  it  at  a 
constant  velocity.  The  wheel  is  moved  longitudinally 
by  displacing  the  axle  node.  This  simulates  a  towed 
wheel  and  duplicates  the  procedure  of  a  rolling  resis¬ 
tance  test  using  the  CRREL  Instrumented  Vehicle. 
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These  models  were  configured  for  both  implicit 
dynamic  analysis  and  explicit  dynamic  analysis.  The 
explicit  models  used  ALE  meshing,  where  the  mesh 
automatically  adjusts  during  the  deformation  by  al¬ 
lowing  mass  transfer  between  elements.  The  final 


Figure  56.  Model  of  a  wheel  rolling  on  20  cm  of 
fresh  snow  (infinite  elements  not  shown). 


mesh,  therefore,  does  not  directly  reflect  the  deforma¬ 
tion  of  each  element,  as  in  a  standard  finite  element 
mesh,  because  the  nodes  and  elements  have  moved 
relative  to  the  snow  mass.  However,  the  final  surface 
shape  is  correctly  represented,  and  tracer  nodes  can 
be  used  to  track  the  mass  movement.  The  ALE  mesh¬ 
ing  significantly  reduces  model  run  time  (from  over 
48  hours  to  approximately  24  hours  on  an  SGI  Ori¬ 
gin),  eliminates  severe  mesh  distortion,  and  smooths 
the  resultant  load  response,  but  it  does  not  affect  the 
numerical  results.  Eigure  57  shows  the  deformed 
mesh  for  the  deep  snow  model  with  and  without  the 
ALE  meshing. 

The  snow  material  properties  were  adjusted  to 
match  field  measurements  of  sinkage  and  motion 
resistance.  This  was  done  primarily  through  changes 
to  the  cap  hardening  and  consequent  adjustments 
to  the  Drucker-Prager  parameters.  Since  snow  com¬ 
paction  is  the  primary  mechanism  contributing  to 
motion  resistance,  the  pressure-volume  relationship  for 


Figure  57.  Comparison  of  standard  mesh  (top)  and  ALE  mesh  (bot¬ 
tom).  (From  Shoop  et  al.  1999.) 
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snow,  reflected  in  the  hardening  table,  has  the  great¬ 
est  influence  on  resulting  tire  forces  and  sinkage. 
Changes  to  the  acceleration,  velocity,  and  friction 
coefficient  have  little  effect  on  model  results,  partly 
because  of  the  low  velocities  and  accelerations  in¬ 
volved  in  standard  snow  rolling  resistance  tests. 

Motion  resistance  forces  and  sinkage 

The  tire-snow  model  was  compared  to  field 
measurement  of  tire  forces  made  with  an  instru¬ 
mented  vehicle.  Model  results  were  also  compared  to 
rolling  resistance  and  sinkage  predictions  made  from 
established  algorithms  for  snow  (as  adopted  by  the 
NATO  Reference  Mobility  Model). 

Experimental  measurements.  Vehicle  performance 
was  measured  using  the  CRREL  Instrumented  Vehi¬ 
cle  (CIV),  which  is  instrumented  to  measure  vertical, 
longitudinal,  and  lateral  forces  at  the  tire-terrain  in¬ 
terface;  wheel  speed  at  each  wheel;  true  vehicle 
speed;  and  steering  angles.  The  CIV  is  designed  as  a 
research  tool  to  perform  various  mobility  tests  (trac¬ 
tion,  resistance,  and  maneuverability)  using  different 
tires,  traction  aids,  and  vehicle  configurations  (of 
braking  and  driving  wheels)  on  a  range  of  terrain 
surfaces  including  dry,  wet,  snow,  ice,  and  freezing 
or  thawing  ground.  The  data  obtained  from  these  tests 
are  used  for  model  validation  and  for  developing  al¬ 
gorithms  to  predict  vehicle  performance  on  cold- 
weather  terrain. 

The  vehicle,  originally  a  1977  American  Motors 
Corporation  Jeep  Cherokee,  is  shown  in  Figure  58. 
Each  wheel  is  instrumented  with  a  three-component 
load  cell  (Fig.  59),  designed  and  calibrated  to  respond 
to  the  forces  at  the  terrain  interface.  The  wheel  speed 
is  measured  using  a  proximity  sensor  with  a  set  of 
100  measurement  pegs  set  into  the  brake  rotor.  Wheel 
speed  is  compared  to  the  true  vehicle  speed  measured 
using  a  fifth  wheel  or  a  sonic  speed  sensor  mounted 


Figure  58.  The  CRREL  Instrumented  Vehicle  (CIV). 
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Figure  59.  Configuration  of  speed 
sensors  and  axle-mounted  load 
cells  on  the  CIV  (From  Shoop  et  al. 

1994.) 

on  the  vehicle.  The  difference  between  the  wheel 
speed  Vw  and  the  vehicle  speed  Vy  is  termed  the  dif¬ 
ferential  interface  velocity  (DIV),  also  called  the 
“slip  speed”  or  “longitudinal  slip  velocity”  (SAE 
1992).  This  is  often  normalized  by  the  speed  of  the 
vehicle  and  called  slip  i  (ISTVS  1977): 

■  =  ^w  -^v  xioo%  .  (28) 

Using  this  convention,  driving  slip  is  positive  and 
braking  slip  is  negative. 

Rolling  resistance  or  motion  resistance  is  the  sum 
of  the  forces  resisting  vehicle  motion.  These  are 
forces  due  to  the  internal  friction  of  the  moving  parts 
of  the  vehicle,  the  air  drag,  the  internal  resistance  of 
the  tires  due  to  flexing  of  the  belts  and  plies  and  vis¬ 
cous  properties  of  the  rubber  compounds,  and  the 
added  resistance  due  to  the  deformation  of  the  sur¬ 
face.  The  CIV  measures  motion  resistance  forces 
directly  at  the  wheel,  eliminating  the  forces  acting  on 
the  vehicle  body.  Motion  resistance  is  tested  with  the 
CIV  by  rolling  the  vehicle  through  snow  at  a  constant 
speed.  The  vehicle  is  driven  with  the  rear  wheels 
only,  and  the  front  brakes  are  released  so  that  the 
only  forces  on  the  front  wheels  are  due  to  the  tires 
and  snow.  The  internal  resistance  of  the  tire  is  meas¬ 
ured  by  rolling  it  on  a  hard  surface.  This  value  is  sub¬ 
tracted  from  the  snow  motion  resistance  measure- 
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ment  to  obtain  the  resistance  due  only  to  the  deform¬ 
able  snow  (or  other  terrain).  The  experimental  values 
of  rolling  resistance  presented  here  are  only  the  por¬ 
tion  of  the  resistance  due  to  the  deformation  of  the 
terrain  (snow). 

CRREL  has  a  large  database  of  vehicle  measure¬ 
ments  in  snow.  For  comparison  to  model  data,  most 
of  the  field  measurements  were  taken  from  a  study  of 
wheeled  and  tracked  vehicle  performance  in  snow 
that  was  performed  in  Houghton,  Michigan,  in  1988 
and  1989  and  documented  in  Blaisdell  et  al.  (1990) 
and  Green  and  Blaisdell  (1991).  Motion  resistance 
measurements  formed  a  significant  portion  of  this 


study.  Resistance  in  shallow  snow  was  studied  fur¬ 
ther  in  an  experimental  program  in  1993,  also  in 
Houghton,  Michigan,  documented  in  Richmond 
(1995).  The  data  from  these  studies  are  given  in  Ta¬ 
ble  11.  The  majority  of  these  data  were  collected  us¬ 
ing  Michelin  XCH4  tires  except  where  a  different 
radius  is  noted.  This  database  was  used  to  develop  a 
semi-empirical  performance  model  for  mobility  on 
snow,  which  was  incorporated  into  the  NATO  Refer¬ 
ence  Mobility  Model  (NRMM  II)  (Richmond  et  al. 
1995,  Ahlvin  and  Shoop  1995).  Both  the  original  data 
and  the  predictive  algorithms  were  compared  to  the 
tire-snow  finite  element  model. 


Table  11.  Measured  sinkage  and  resistance  in  snow.  All 
tires  are  the  same  (radius  =  38  cm),  except  for  the  last  four 
measurements. 


Snow  depth 
(m) 

Sinkage 

(m) 

Resistance 

coefficient 

Snow  density 
(kg/rrP) 

0.13 

0.10 

0.05 

160 

0.16 

0.12 

0.06 

160 

0.15 

0.12 

0.05 

170 

0.19 

0.15 

0.08 

170 

0.16 

0.13 

0.07 

160 

0.17 

0.14 

0.06 

160 

0.18 

0.15 

0.06 

170 

0.19 

0.16 

0.07 

170 

0.21 

0.17 

0.07 

220 

0.21 

0.17 

0.10 

220 

0.20 

0.16 

0.07 

220 

0.21 

0.17 

0.12 

220 

0.20 

0.16 

0.07 

220 

0.21 

0.16 

0.09 

220 

0.18 

0.14 

0.07 

220 

0.19 

0.16 

0.10 

220 

0.21 

0.16 

0.06 

220 

0.22 

0.17 

0.07 

220 

0.21 

0.15 

0.07 

220 

0.18 

0.13 

0.05 

220 

0.18 

0.13 

0.06 

220 

0.20 

0.15 

0.07 

220 

0.21 

0.16 

0.07 

220 

0.22 

0.17 

0.08 

220 

0.22 

0.17 

0.06 

220 

0.22 

0.17 

0.05 

220 

0.14 

0.10 

0.02 

110 

0.13 

0.90 

0.03 

110 

0.17 

0.14 

0.03 

110 

0.16 

0.12 

0.04 

110 

0.18 

0.15 

0.03 

110 

0.16 

0.12 

0.06 

110 

0.16 

0.13 

0.14 

230 

0.38 

0.33 

0.19 

220 

0.38 

0.31 

0.19 

160 

0.36 

0.31 

0.21 

160 

0.32  (r  =  35  cm) 

0.26 

0.17 

190 

0.28  (r  =  33  cm) 

0.25 

0.10 

120 

0.26  (r  =  33  cm) 

0.22 

0.16 

170 

0.25  (r  =  33  cm) 

0.23 

0.12 

120 
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Figure  60.  Parameters  used  to  predict  motion  resistance  us¬ 
ing  the  NRMM  aigorithm.  (From  Richmond  1995.) 


NRMM  prediction  of  motion  resistance  in  snow. 
The  following  method  was  proposed  by  CRREL  for 
predicting  the  motion  resistance  in  snow.  It  was 
adopted  into  the  NATO  Reference  Mobility  Model  as 
part  of  the  overall  vehicle  performance  predictions 
scheme.  Only  the  portions  of  the  algorithms  related  to 
motion  resistance  in  snow  are  stated  here.  Based  on 
Richmond  (1995),  the  motion  resistance  for  the  lead¬ 
ing  tire  in  shallow  snow  is  based  on  the  vehicle 
sinkage  z  according  to 

=13.6041(poWa)‘'^®  (29) 

a  =  r  arcos  [(r-z)/r]  (30) 


z  = 


1-- 


Po 


0.519  +  0.0023/7^ 


(31) 


where 

Po 

w 

a 

r 

z 


PmdLX 


motion  resistance  (leading  tire  only)  (N) 

initial  snow  density  (kg/m^) 

maximum  tire  width  (m) 

arc  length  in  contact  with  snow  (m) 

tire  radius  (m) 

sinkage  (m) 

snow  depth  (m) 

maximum  tire  contact  pressure  (kPa). 


These  parameters  are  illustrated  in  Figure  60. 

Very  little  data  is  available  for  wheeled  vehicles  in 
deep  snow.  However,  based  on  numerous  observa¬ 
tions  of  snow  vehicle  movement,  the  following  fac¬ 
tors  are  applied  to  the  shallow  snow  equations  to  es¬ 
timate  the  additional  motion  resistance  in  deep  snow. 
The  factors  represent  an  engineering  estimate  of  the 
additional  forces  due  to  plowing  of  the  vertical  face 
of  the  wheel  and  undercarriage  drag  (Richmond  et  al. 
1995).  The  deep  snow  modifiers  would  be  applicable 
to  wheels  being  dragged  or  towed  through  deep 


snow,  since  a  driven  wheel  would  not  likely  have 
enough  traction  to  overcome  the  large  resistance  and 
would  become  immobile. 

For  Po  >  150  kg/w?  and  2/3  r  <  z  <  r,  (32) 

Factor  =1.5 

For  Po  >  150  kg/w?  and  z  >  r,  (33) 

Factor  =  2.5. 

Comparison  of  measurements,  NRMM,  and  finite 
element  results.  In  the  FEA  model  the  longitudinal 
reaction  force  at  the  wheel  hub  represents  the  motion 
resistance  force  on  the  wheel  due  to  deformation  of 
the  snow  and  is  directly  comparable  to  the  motion 
resistance  measured  with  the  CRREF  Instrumented 
Vehicle.  Similarly  the  vertical  displacement  of  the 
hub  node  is  equivalent  to  the  sinkage  of  the  wheel 
into  the  snow  (since  the  wheel  itself  is  not  deforming 
in  this  simulation).  Figure  61  shows  the  resistance 
coefficient  (the  longitudinal  load  divided  by  the  verti¬ 
cal  load  on  the  hub)  and  vertical  displacement  for  one 
of  the  snow  model  runs.  For  this  model  the  vertical 
loading  occurs  from  0  to  1.5  s,  the  horizontal  accelera¬ 
tion  occurs  from  1.5  to  2.0  s,  and  then  constant  speed 
translation  is  maintained  from  2.5  to  5  s  at  a  speed  of  8 
kph  (5  mph).  The  sinkage  and  motion  resistance  val¬ 
ues  are  chosen  as  the  average  of  the  values  occurring 
during  the  constant  speed  portion  of  the  simulation. 

Comparisons  of  the  measured  and  modeled  sink- 
age  and  motion  resistance  are  shown  in  Figure  62  as 
a  function  of  snow  depth.  The  model  overestimates 
motion  resistance,  but  the  trends  are  consistent  with 
the  measured  data  and  the  model  falls  directly  along 
the  upper  edge  of  the  measurements,  constituting  a 
conservative  prediction.  The  NRMM  predictions  of 
motion  resistance  are  discontinuous  based  on  the 
snow  depth.  For  this  tire,  motion  resistance  for  snow 
depths  less  than  25  cm  are  calculated  using  NRMM 
eq  29,  for  depths  from  25  to  38  cm  using  eq  31,  and  for 
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Figure  61.  Modeled  resistance  force  in  coefficient  form  (longitudi¬ 
nal/vertical)  and  vertical  displacement  of  the  wheel  hub  for  the  FEM 
wheel  rolling  through  20-cm  snow. 


Figure  62.  Finite  element  model,  measured  data,  and  NRMM 
predictions  for  sinkage  and  motion  resistance  in  fresh  snow. 
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snow  depths  greater  than  38  cm  using  eq  32.  The 
finite  element  model  falls  between  the  NRMM  mo¬ 
tion  resistance  prediction  equations  for  deep  snow,  as 
does  the  measured  data.  The  bottom  part  of  Figure  62 
shows  the  corresponding  sinkage,  which  the  finite 
element  model  (with  zero  slip)  underestimates, 
particularly  for  the  deeper  snow. 

The  NRMM  predictions  are  also  a  function  of  the 
initial  snow  density.  Although  the  material  model 
was  designed  to  simulate  200-kg/m^  snow,  NRMM 
predictions  based  on  240-kg/m^  snow  yield  an  excel¬ 
lent  match  to  the  finite  element  resistance,  as  seen  in 
the  top  of  Figure  63.  This  would  indicate  that  the 
FEA  model  is  more  representative  of  a  slightly 
denser  snow. 

Figure  63  also  shows  model  results  for  a  wheel 


with  zero  slip  and  with  unrestricted  slip.  The  NRMM 
deep  snow  resistance  predictions  are  close  to  the  fi¬ 
nite  element  results  for  a  free  slipping  wheel,  as 
would  be  expected  for  a  wheel  being  dragged  through 
deep  snow,  while  shallow  snow  results  are  best  repre¬ 
sented  by  the  simulation  with  zero  slip. 

The  bottom  of  Figure  63  shows  the  measured 
sinkage  data  falling  between  the  zero  slip  and  unre¬ 
stricted  slip  models.  Measured  slip  in  shallow  snow  is 
near  zero,  but  higher  slip  measurements  are  not  un¬ 
common  as  the  snow  gets  deeper.  Realistically  the 
model  should  behave  somewhere  between  these  two 
extremes.  This  part  of  the  physics  of  the  interaction  is 
not  accurately  reproduced  in  the  model,  partly  be¬ 
cause  the  tread,  which  is  designed  to  resist  slippage, 
is  not  yet  accurately  modeled. 


Figure  63.  Finite  eiement  simuiations  at  zero  siip  and  at  uniimited 
siip  and  NRMM  motion  resistance  predictions  for  240-kg/m^  snow. 
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Snow  deformation  beneath  the  wheel 

Snow  deformation  under  a  vehicle  can  be  meas¬ 
ured  by  excavating  a  cross  section  of  the  snow  as 
shown  earlier  in  Figure  3.  Additional  detail  can  be 
obtained  by  marking  the  snow  prior  to  deformation 
by  backfilling  a  small-diameter,  vertical  hole  with 
dark-colored  chalk  dust.  The  vehicle  is  then  driven 
into  the  snow,  and  the  changes  in  the  lines  can  be 
seen  when  the  cross  section  is  excavated.  This  proc¬ 
ess  is  illustrated  in  Figure  64. 

Once  the  snow  cross  section  has  been  excavated, 
snow  density  can  be  measured  by  collecting  and 
weighing  samples  of  a  known  volume.  This  was  done 
in  several  experiments  (Richmond  1995).  Modeled 
deformation  of  20-cm,  200-kg/m^  snow  was  com¬ 
pared  to  measured  snow  deformation  beneath  the 
CIV  in  similar  snow  conditions.  Because  the  model 
used  ALE  meshing,  the  displacement  of  the  snow 
was  documented  by  placing  tracer  particles  in  the 
model,  spaced  at  the  same  distance  as  the  markings  in 
the  field.  Comparisons  of  the  measured  and  modeled 


snow  deformation  in  cross  sections  in  the  direction  of 
travel  (longitudinal)  and  transverse  to  the  direction  of 
travel  are  shown  in  Figures  65  and  66.  The  deviatoric 
stress  contours  shown  in  the  longitudinal  cross  sec¬ 
tion  reveal  the  location  where  a  new  shearing  surface 
begins  to  develop  ahead  of  the  advancing  wheel.  The 
failure  arc,  even  for  snow,  clearly  advances  down¬ 
ward  and  away  from  the  wheel,  just  as  shown  ex¬ 
perimentally  for  sand  in  Figure  1.  The  deformation 
predicted  in  the  model  is  slightly  less  than  the  ob¬ 
served  deformation,  partly  because  the  model  is  for 
deeper  (20  cm)  and  denser  (200  kg/m^)  snow  than 
occurred  in  the  field.  The  general  shape  of  the  dis¬ 
placement  is  the  same,  however. 

The  measured  density  and  the  modeled  density  are 
in  good  agreement,  as  shown  in  the  transverse  cross 
section  in  Figure  67.  The  modeled  density  directly 
beneath  the  wheel  (460  kg/m^)  is  only  10%  different 
from  the  measured  density  (510  kg/m^).  Near  the  tread 
shoulder,  the  modeled  density  was  290-380  kg/m^, 
with  measured  values  ranging  from  280  to  300  kg/m^. 


Figure  64.  Marking  the  snow  to  observe  snow  de¬ 
formation  after  vehicie  passage. 
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Viewport:  1  Model:  Model-1  Module:  Visualization 


Scale  Factor;  +1.00 


S,  Mises 
(Ave.  Crit.:  75%) 


Figure  65.  Comparison  of  measured  dispiacement  (14.5-cm  snow)  to 
modeied  dispiacement  and  deviatoric  stress  (20-cm  snow)  in  the  ion- 
gitudinai  direction. 


Figure  66.  Comparison  of  measured  (19-cm  snow)  and  modeied  (20-cm  snow)  dispiacement  in  a  cross 
section  transverse  to  the  direction  of  travei. 
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RESTART  FILE  =  20cm2-3 
TOTAL  TIME  3.07 


Figure  67.  Comparison  of  measured  (19-cm  snow)  and 
modeied  (20-cm  snow)  densities  in  a  cross  section  trans¬ 
verse  to  the  direction  of  travei. 


Deformable  tire  on  soil 

On  harder  terrain  it  is  important  to  model  the  de¬ 
formation  of  the  tire  in  addition  to  the  deformation  of 
the  terrain.  Because  of  the  large  deformation  and  run 
times  involved,  it  is  important  that  the  tire  model  be 
as  efficient  as  possible  while  still  maintaining  an  ac¬ 
curate  contact  patch.  Thus,  the  tire  model  used  for 
rolling  on  soil  was  based  on  the  formulation  proposed 
by  Darnell  et  al.  (1997). 

The  Shoop-Darnell  tire  model  was  modified  to 
roll  across  deformable  material,  such  as  soil.  This  tire 
was  first  placed  on  an  elastic  terrain  material  with  the 
properties  of  compacted  sand.  The  simulation  was 
conducted  in  four  steps:  the  tire  was  first  inflated, 
then  lowered  into  the  soil,  allowed  to  come  to  equi¬ 
librium,  and  then  rolled  to  the  left,  as  seen  in  Figure 
68.  Only  the  tread  and  soil  elements  are  shown  in  the 
figure;  the  user-defined  sidewall  elements  are  not 
displayed.  Figure  68d  is  an  oblique  view  of  the  three- 
dimensional  model  showing  the  contact  pressure. 


The  same  tire  was  then  rolled  on  a  Drucker-Prager 
cap  plasticity  model  representing  the  McCormick 
Ranch  sand  (discussed  earlier).  The  results  of  this 
simulation  are  given  in  Figure  69,  showing  the  maxi¬ 
mum  principal  plastic  strain  and  the  principal  stress 
directions.  The  maximum  principal  stress  is  compres¬ 
sive  beneath  and  slightly  forward  of  the  tire.  The 
minimum  principal  stress  is  tensile  and  is  largest  just 
below  and  to  the  sides  of  the  tire,  oriented  away  from 
the  tire.  The  difference  between  these  stresses  indi¬ 
cates  the  magnitude  of  shear  stress  in  the  soil,  which 
is  largest  beneath  and  directly  in  front  of  the  tire.  This 
corresponds  to  the  general  shear  zones  beneath  a 
towed  rigid  wheel  observed  by  Wong  and  Reece 
(1967),  shown  in  Figure  1. 

The  tire-soil  models  are  prototypes  requiring  re¬ 
finement  and  development.  They  will  eventually  be 
compared  to  field  measurements  of  tire  forces,  de¬ 
flection,  and  contact  area,  including  sinkage  and 
deformation  of  the  soil  surface. 
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a.  Tire  inflated. 


b.  Tire  lowered  onto  soil. 


c.  Tire  rolled  to  the  left. 

Figure  68.  Rolling  tire  on  an  elastic  material  (sand). 
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d.  Contact  pressure. 
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Figure  68  (cont.). 


Viewports  Model:  Model-1  Module:  Visualization 
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a.  Plastic  strain. 


b.  Principal  stress  directions  as  the  tire  rolls  to  the  left. 


Figure  69.  Rolling  Shoop-Darnell  tire  on  a  Drucker-Prager 
cap  model  of  the  McCormick  Ranch  sand. 
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Figure  70.  Deformed  mesh  of  the  modal  analysis  tire  model  in 
snow  after  8  cm  of  sinkage  in  snow. 


Figure  71.  Close-up  of  a  tire  after  sinking  8  cm  into  the 
snow  at  the  beginning  of  roll.  The  tire  shows  very  little  de¬ 
formation  (the  deflection  is  3  mm),  indicating  that  a  rigid 
wheel  may  be  a  good  approximation. 


Deformable  tire  on  snow 

Simulations  of  the  modal  analysis  tire  on  snow  are 
in  progress.*  Although  preliminary  runs  are  very 
slow  and  have  not  progressed  to  completion,  inter¬ 
mediate  results  are  available.  Figure  70  shows  the 
simulation  after  the  tire  has  been  lowered  onto  the 
snow  and  begins  to  roll.  The  sinkage  is  approxi¬ 
mately  8  cm,  which  is  about  the  same  as  the  modeled 
sinkage  for  the  rigid  wheel  on  20-cm  snow,  as  seen  in 
Figure  61.  Of  significant  note,  however,  is  that  the 
deformable  tire  has  undergone  very  little  deformation 
(the  deflection  is  3  mm,  which  is  less  than  2%),  as 


*  Personal  communication  with  K.  Kestler,  private 
consultant,  2000. 


seen  in  Figure  71,  revealing  that  the  rigid  wheel 
model  for  deep  fresh  snow  may  be  a  valid  assumption 
(although  these  results  are  only  preliminary). 

CONCLUSIONS 

Summary 

The  interaction  between  a  tire  and  deformable  ter¬ 
rain  is  a  complicated  dynamic  process  that  involves 
the  deformation  of  both  the  tire  and  the  terrain  in 
three  dimensions.  The  tire  is  a  complex  structure 
made  of  many  materials,  and  the  terrain  undergoes 
large  deformations  in  an  inelastic  manner,  so  neither 
component  is  easily  modeled.  The  objective  of  this 
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research  was  to  produce  a  three-dimensional  finite 
element  model  of  tire-terrain  interaction  that  can  be 
used  to  explore  the  effects  of  tire  and  terrain  variables 
on  vehicle  mobility.  Such  a  model  can  be  used  for 
tire  design  and  specification  for  off-road  vehicles,  for 
vehicle  performance  prediction,  and  for  terrain  dam¬ 
age  prediction  and  reduction  of  the  environmental 
impact  of  off-road  travel. 

The  details  of  the  tire-soil  modeling  problem  are 
divided  into  three  topic  areas:  1)  material  models  for 
the  terrain  material,  2)  tire  models  for  use  on  a  de¬ 
formable  substrate,  and  3)  the  combined  tire-terrain 
model  and  the  treatment  of  the  interface.  The  terrain 
materials  simulated  in  this  study  were  fresh  snow  and 
compacted  sandy  soil.  The  two  material  models  used 
were  a  modified  Drucker-Prager  cap  plasticity  model 
and  a  critical-state,  crushable  foam  model.  Both 
models  were  considered  suitable  for  capturing  the 
highly  compressible  behavior  of  fresh,  natural  snow 
(initial  density  of  200-250  kg/m^  at  temperatures  of 
-10°  to  -1°C).  Model  parameters  were  generated 
from  field  test  data  and  from  the  literature,  matching 
this  snow  type  as  closely  as  possible.  The  snow 
model  was  validated  using  plate  sinkage  test  data  for 
snow  of  similar  age  and  density.  The  soil  model 
represents  sand  similar  to  that  used  during  vehicle 
mobility  experiments  at  CRREL.  The  material  was 
modeled  using  a  Drucker-Prager  cap  plasticity  model 
with  input  parameters  from  the  literature. 

To  apply  a  tire  model  to  deformable  terrain,  the 
model  must  be  efficient  yet  accurately  portray  the  tire 
structural  behavior.  Four  tire  models  were  evaluated 
for  suitability  to  rolling  on  deformable  terrain:  1)  a 
rigid  tire  model,  2)  a  simplified  tire  model  using 
methodology  developed  by  Darnell  at  the  University 
of  Michigan  for  use  in  vehicle  dynamics  simulations, 
3)  a  tire  model  of  the  type  used  for  harmonic  vibra¬ 
tion  modal  analysis,  with  a  smooth  tread,  and  4)  a  tire 
model  similar  to  model  3  except  with  a  straight 
ribbed  (longitudinal)  tread.  All  of  the  models  were 
built  to  represent  tires  used  in  the  experimental  test 
program  for  comparison  to  measured  tire  behavior  in 
terms  of  deflection,  contact  area,  deflected  sidewall 
profile,  contact  stress  distribution,  and  rolling  resis¬ 
tance  forces  on  deformable  terrain  (snow  and  soil). 

Three  models  of  the  combined  tire  and  terrain 
were  developed.  The  first  is  a  rigid  tire  on  fresh 
snow.  The  second  is  the  Shoop-Darnell  tire  model 
rolling  on  a  soil,  and  the  last  is  the  modal  analysis  tire 
on  snow.  Model  results  simulating  a  rigid  tire  rolling 
on  snow  were  compared  to  tire  forces  measured  using 
an  instrumented  vehicle.  The  measured  snow  defor¬ 
mation  under  the  wheel  was  also  compared  to  model 
results.  The  model  exhibits  good  agreement  with  mo¬ 


tion  resistance  forces  but  underestimates  sinkage.  The 
measured  deformation  patterns  are  duplicated  in  the 
model,  but  the  modeled  displacements  are  less  than 
what  is  measured  in  the  field.  Model  results  were  also 
compared  to  snow  mobility  predictions  made  using 
the  largely  empirical  NATO  Reference  Mobility 
Model,  with  good  agreement  for  forces  but  underes¬ 
timation  of  sinkage.  The  amount  of  wheel  slip  had  a 
major  impact  on  the  modeled  sinkage. 

Simulations  of  a  tire  rolling  on  soil  and  snow  us¬ 
ing  the  Shoop-Darnell  tire  and  the  modal  analysis 
tire  model  are  operational  but  have  not  been  vali¬ 
dated.  Of  significant  note,  however,  is  that  a  simula¬ 
tion  using  the  modal  analysis  tire  on  fresh  snow  pre¬ 
dicts  very  little  tire  deformation  (3  mm,  or  less  than 
2%  deflection).  This  suggests  that  the  assumption  of 
a  rigid  tire  may  be  used  for  soft  terrain  such  as  deep, 
fresh  snow  without  loss  of  model  accuracy. 

Significant  findings 

The  following  summarizes  the  significant 
achievements  of  this  work: 

1.  A  material  model  was  developed  for  fresh  snow 
and  validated  with  plate  sinkage  tests  in  the  lab 
and  field.  Good  agreement  of  measured  and  mod¬ 
eled  forces,  displacements,  and  changes  in  mate¬ 
rial  density  were  achieved. 

2.  Evaluations  of  several  finite  element  tire  models 
suitable  for  rolling  on  a  deformable  substrate  in¬ 
dicate  that  the  Darnell  model  yielded  accurate  re¬ 
sults  and  was  computationally  efficient.  The  mo¬ 
dal  analysis  type  of  tire  model  is  also  suitable  but 
is  more  computationally  intensive. 

3.  Combined  tire-terrain  models  utilizing  the 
Shoop-Darnell  tire  and  the  modal  analysis  tire 
model  are  operational  but  have  not  yet  been  vali¬ 
dated.  Preliminary  results  of  the  modal  analysis 
tire  model  on  snow  show  very  little  deformation 
in  the  tire,  indicating  that  the  rigid  wheel  simplifi¬ 
cation  may  be  valid  for  soft  terrain. 

4.  A  model  of  a  rigid  wheel  on  fresh  snow,  validated 
experimentally,  shows  good  agreement  with 
measured  motion  resistance  forces,  snow  dis¬ 
placement,  and  snow  compaction  and  agrees  with 
results  predicted  by  the  NATO  Reference  Mobil¬ 
ity  Model. 

5.  The  rigid  wheel  on  snow  model  does  not  capture 
the  tread  effects  of  the  tire-terrain  interaction  (as 
tread  patterns  were  not  expressly  modeled).  For  a 
free-rolling  wheel,  the  amount  of  slip  and  sinkage 
is  overpredicted.  When  the  model  wheel  slip  is  set 
to  zero,  as  is  measured  during  shallow  snow  resis¬ 
tance  tests  in  the  field,  sinkage  is  underpredicted. 
The  most  realistic  results  lie  between  these  two 
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cases  and  can  possibly  be  modeled  with  “rough” 
friction  attributes.  A  full  understanding  of  this  ef¬ 
fect  will  be  critical  for  traction  and  braking  stud¬ 
ies. 

Recommended  applications  and  future 
research  needs 

Parametric  design  anaiysis 

The  primary  applications  of  a  tire-terrain  model 
are  1)  all-season  and  off-road  performance  prediction 
for  tires,  2)  tire  design,  optimizing  geometry,  materi¬ 
als,  and  tread  patterns,  and  3)  specification  of  the 
correct  tires  based  on  expected  operating  conditions. 
Naturally,  design  use  of  the  model  could  include  a 
parametric  analysis  of  the  impact  of  tire  and  terrain 
parameters  on  tire  performance.  Such  a  study,  evalu¬ 
ating  the  effects  of  tire  geometry  and  loading  on  roll¬ 
ing  resistance  in  snow,  may  now  be  undertaken  using 
the  results  of  this  research. 

Interface  friction  and  interlocking  of 
tread  and  terrain 

The  processes  occurring  at  the  interface  of  the  tire 
and  terrain  are  also  of  great  interest.  Advances  in 
driving  and  braking  traction  mechanics  through  so¬ 
phisticated  study  of  interface  friction,  and  the  imple¬ 
mentation  of  this  in  a  numerical  model,  are  of  major 
significance  because  of  their  large  impact  on  tire  and 
vehicle  performance  and  safety.  Since  the  tread  pat¬ 
tern  was  not  modeled  in  this  study,  the  details  of  the 
interface  friction  and  frictional  interlocking,  which 
affected  the  resulting  slip  of  the  tire,  were  not  spe¬ 
cifically  modeled.  Issues  include  the  interaction  of 
terrain  material  and  tread  blocks  and  the  impacts  of  a 
spinning  tire  on  traction  through  frictional  heating 
and  mass  shearing.  Modeling  challenges  include  cou¬ 
pling  with  other  models,  such  as  submodels  of  tread 
blocks,  and  discrete  element  or  hydrodynamic  models 
of  terrain  material.  The  practical  application  of  this  is 
readily  apparent  for  tread  design. 

Deformation  of  terrain  materials 
subjected  to  vehicle  loading 

Research  is  also  needed  in  areas  relating  to  the  de¬ 
formation  of  the  materials  beneath  the  wheel,  such  as: 

1.  Soil  deformation  under  vehicles  and  implications 
for  terrain  damage  through  soil  compaction,  rut 
formation,  root  tearing,  etc. 

2.  Washboard  formation  on  secondary  roads,  includ¬ 
ing  the  effects  of  suspension  parameters,  vehicle 
speed,  and  road  material. 

3.  Dynamic  models  of  vehicles  on  pavements,  in¬ 
cluding  impacts  on  pavement  structure  and  dete¬ 


rioration,  loading  effects  on  soil  water  (pumping 
up  through  pavement  layers  and  cracks),  and  the 
function  of  drainage  materials,  including  geosyn¬ 
thetics. 

4.  The  impact  of  temperature-dependent  material 
properties,  such  as  frictional  heating  and  melting 
of  snow,  ice,  or  frozen  ground. 

Additional  terrain  material  models 

Future  work  should  also  include  extending  the  ma¬ 
terial  models  to  other  terrain  materials,  such  as  wet, 
trafficked  snow  or  slush,  and  soft,  loose  soils.  Traf¬ 
ficked  snow  is  typical  of  much  of  the  snow  on  roads 
before  clearing  and  therefore  is  of  primary  interest  to 
snow  maintenance  operations,  tire  companies  for 
improving  design  of  snow  tires,  automotive  compa¬ 
nies  for  optimizing  vehicle  handling,  and  for  the  roll¬ 
ing  resistance  of  the  snow  and  its  effect  on  fuel  econ¬ 
omy.  The  thawing  soils  material  represents  a  critical 
case  for  roads  subjected  to  freeze-thaw,  where  soil 
moisture  is  pulled  to  the  surface  during  freezing  (ex¬ 
panding  soil  pores  or  forming  ice  lenses),  resulting  in 
a  wet  and  loose  soil  during  seasonal,  intermittent,  or 
even  daily  thawing  events.  Some  of  the  difficulties  in 
modeling  this  unique  material  may  be  simulated  us¬ 
ing  the  more  flexible  Multi-Mechanical  Model  pro¬ 
posed  by  Peters*  and  Smith  (2000)  or  through  hydro- 
dynamic  modeling.  Additional  behaviors  specific  to 
the  interaction  of  snow  grains  during  loading  and 
resulting  sintering  may  be  better  approached  using  a 
discrete  element  model  as  proposed  by  Johnson  and 
Hopkins.** 
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